Identification of novel tumor microenvironment-associated genes in gastric cancer based on single-cell RNA-sequencing datasets

Abstract

Tumor microenvironment and heterogeneity play vital roles in the development and progression of gastric cancer (GC). In the past decade, a considerable amount of single-cell RNA-sequencing (scRNA-seq) studies have been published in the fields of oncology and immunology, which improve our knowledge of the GC immune microenvironment. However, much uncertainty still exists about the relationship between the macroscopic and microscopic data in transcriptomics. In the current study, we made full use of scRNA-seq data from the Gene Expression Omnibus database (GSE134520) to identify 25 cell subsets, including 11 microenvironment-related cell types. The MIF signaling pathway network was obtained upon analysis of receptor–ligand pairs and cell–cell interactions. By comparing the gene expression in a wide variety of cells between intestinal metaplasia and early gastric cancer, we identified 64 differentially expressed genes annotated as immune response and cellular communication. Subsequently, we screened these genes for prognostic clinical value based on the patients’ follow-up data from The Cancer Genome Atlas. TMPRSS15, VIM, APOA1, and RNASE1 were then selected for the construction of LASSO risk scores, and a nomogram model incorporating another five clinical risk factors was successfully created. The effectiveness of least absolute shrinkage and selection operator risk scores was validated using gene set enrichment analysis and levels of immune cell infiltration. These findings will drive the development of prognostic evaluations affected by the immune tumor microenvironment in GC.

Introduction

Although the incidence and mortality rates have been declining worldwide, gastric cancer (GC) remains a common and lethal malignancy, especially in Asian countries (Smyth et al., 2020). In addition to traditional chemotherapy and surgery, adjuvant therapies, such as molecularly targeted therapy (Mundekkad and Cho, 2022) and immunotherapy, are gradually emerging as the staples of GC treatment. Uncovering the molecular mechanisms of the initiation and progression of GC is therefore critical for improving therapeutic efficacy. Molecular markers help deepen our understanding of GC subtypes (Bijlsma et al., 2017; Zhao et al., 2019), and with the technological advances in high-throughput sequencing, the focus on tumor heterogeneity is driving progress in precision medicine (Zeng and Jin, 2021) simultaneously. There is growing evidence that intra-tumoral heterogeneity includes not only genomic features but also the complex tumor microenvironment (TME). It is well recognized that TME comprises various stromal cells, abundant angiogenesis, and immune cell infiltration (Hanahan, 2022). Understanding the “soil” on which the “seed” grows into a tumor has essential implications for the diagnosis and treatment of GC.

Single-cell RNA sequencing (scRNA-seq) is a potent tool to obtain complete RNA transcripts at the level of single cells by RNA extraction, reverse transcription, amplification, and sequencing (Tang et al., 2009; Ramskold et al., 2012; Picelli et al., 2014). Compared with traditional sequencing methods of tumor tissue, scRNA-seq solves the problem that individual differences between cell types are ignored when the expression of all the genes in thousands of cells is averaged. It shows clear superiority in studying the diversity of tumor cell lineages and predicting interactions between cancer and the microenvironment (Muller and Diaz, 2017; Kumar et al., 2021). In terms of GC research, scRNA-seq shed light on the transcriptome network at different stages of the disease process, from atrophic gastritis, intestinal metaplasia (IM), and dysplasia, to early gastric cancer (EGC) (Zhang et al., 2019), as well as the spatial heterogeneity of microenvironment-related cells in diffuse-type GC (Jeong et al., 2021). Performing scRNA-seq of metastatic GC, the origins of transcriptomic heterogeneity in peritoneal carcinomatosis were analyzed (Wang R. et al., 2021), identifying CLDN11 and CDK12 as markers of lymph node metastasis (Wang B. et al., 2021). Nonetheless, combining large amounts of scRNA-seq data with multi-omics datasets and determining the clinical implications remains a challenge (Wang et al., 2022).

We data-mined existing single-cell transcriptome data from the Gene Expression Omnibus (GEO) to identify microenvironment-related cell types and draw signaling pathway networks according to the marker genes or intercellular communication-related genes. Enrichment analysis was performed to reveal the function of differentially expressed genes (DEGs) between IM and EGC. We leveraged clinical data from The Cancer Genome Atlas (TCGA) to establish a risk-scoring model by the LASSO-Cox regression algorithm. It contained the expression levels of TMPRSS15, VIM, APOA1, and RNASE1, which were highly correlated with clinical outcomes and immune cell infiltration. The results of our study may provide new options for prognostic biomarkers in GC.

Materials and methods

Data collection and pre-processing

ScRNA-seq data of normal cells and GC cells from the GSE134520 dataset were downloaded from the GEO database from the official website (https://www.ncbi.nlm.nih.gov/geo/) (Zhang et al., 2019), which had been pre-processed using CellRanger software. The dataset contains Homo sapiens samples sequenced by the HiSeq X Ten platform. A total of 13 biopsy specimens from the gastric antral mucosa of patients with non-atrophic gastritis, chronic atrophic gastritis, IM, and EGC were selected, and 56,440 cells were included in the dataset. Cells with an abnormal proportion of mitochondrial genes should be removed because it reflects the imbalance of cellular homeostasis and low cell quality. Given the potential existence of diploid cells, cells with genes < 200 or >5,800 were filtered out. Furthermore, the scDblFinder function in the Seurat R package was called to remove the double droplets. Finally, 48,566 cells were enrolled in our study.

The R package “TCGAbiolinks” (version 2.22.2) was used to obtain the gene expression in units of fragments per kilobase million (FPKM) for 407 patients with stomach adenocarcinoma (STAD) from TCGA. We also downloaded the clinical data for these patients, which include survival status (367 patients), stages, age, grades, and survival time. At the same time, the data of “Masked Somatic Mutation” and “Masked Copy Number Segment” called MuTect2 were downloaded.

ScRNA-seq analysis using Seurat

First, we installed the R packages in R (version 4.1) and Seurat (version 4.0.5) (Satija et al., 2015), and used the merge function to merge the created Seurat object. To reduce the influence of different sequencing depths in cells, the normalization for raw counts was performed using the NormalizeData function (“LogNormalize” method), the top 2,000 variable features were identified using the Find Variable Genes function (“vst” method), and the data were integrated and scaled using ScaleData function. Subsequently, we performed principal component analysis (PCA) with the variable genes as inputs and identified the significant principal components whose p-value distributions were then visualized by the jackStraw function. In the FindNeighbors and FindClusters function with a resolution of 0.8, the Louvain algorithm was chosen for cell clustering. FindAllMarkers function with the Wilcoxon rank-sum test was carried out to identify specific marker genes, which compared expression values between cells in the cluster and all other cells. Finally, the results are represented with tSNE (t-distributed Stochastic Neighbor Embedding) dimension reduction by RunTSNE function.

Cell annotation and identification of differentially expressed genes

Cell types were identified based on marker gene sets listed in Table 1. Subsequently, we used the FindAllMarkers function in Seurat to detect DEGs among the different cell types, using a p-value < 0.05 and |log2FoldChange (FC) | > 1 as the thresholds. Additionally, we used the FindAllMarkers function to analyze the differences between the EGC and IM groups to extract DEGs using the same thresholds previously stated. The intersection of three methods (“wilcox,” “t,” and “roc” test) was taken and transformed into a list of the final DEGs. Finally, the expression patterns of these DEGs in the different cell types were shown using heat maps, and the corrplot R package (version 0.92) was then used to perform a correlation analysis between these DEGs.

Enrichment analysis of differentially expressed genes between early gastric cancer and IM groups

Gene Ontology (GO) enrichment analysis is commonly performed in large-scale functional enrichment analyses of genes on different dimensions and levels. It is generally performed on three levels, namely, biological process (BP), molecular function (MF), and cellular component (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database dealing with biological systems such as genomes, biological pathways, diseases, and drugs. The clusterProfiler (version 4.2.0) R package was used to perform GO functional annotation and KEGG pathway enrichment analysis for the DEGs between the EGC and IM groups to identify significantly enriched biological processes. The enrichment results were visualized using bar graphs and bubble graphs. An adjusted p-value of < 0.05 was defined as the significance threshold in the enrichment analysis.

Analysis of cell-cell communication

The CellChat (http://www.cellchat.org/) R package (Jin et al., 2021) was used to calculate the intensity of cell–cell interactions and communication based on single-cell gene expression profiles and known ligands, receptors, and their cofactors. Significant ligand-receptor pairs were further identified based on the probability of receptor–ligand interactions and the results of perturbation testing. We then built a cell–cell communication network by adding the number or the intensity of significantly interacting ligand–receptor pairs between cell types. In addition, we visualized multiple ligand–receptor pairs or intercellular communications mediated by signaling pathways using bubble graphs to investigate the intensity of ligand–receptor interactions between cell types or the characteristics of ligands and receptors in terms of gene expression levels, including their commonalities or differences. Finally, we performed a systemic analysis of the cell–cell communication network by identifying the pathways that contributed most to the incoming and outgoing signals for specific cell groups, performing a network centrality analysis of the identified pathways, and then visualizing the network centrality scores.

Mutation analysis

With the help of GenePattern (https://cloud.genepattern.org), the data of somatic mutations and copy number variations (CNVs) downloaded from TCGA database were analyzed by GISTIC 2.0 to assess the CNV events at the chromosomal arm level and the minimum common region between samples. Then, we use the maftools R package (Mayakonda et al., 2018) to visualize the aforementioned analysis results in the mutation annotation format, and the plotmafSummary function was used to plot the summary file.

Construction of the prognostic model and clinical statistical analysis for model evaluation

To identify genes associated with the prognosis of gastric cancer, we selected the DEGs between the EGC and IM groups as candidate genes and performed univariate Cox regression analysis (p value < 0.01). Forest plots were created (R package: ggplot2) to display each variable’s p-value, hazard ratio (HR), and 95% confidence interval. We then performed least absolute shrinkage and selection operator (LASSO) regression analysis based on these prognostic genes to construct a prognostic model (R package: Glmnet, survival). The risk score was calculated by the formula:

Here, the coef (k) represents the LASSO-Cox regression coefficient, n represents the number of genes, and x(k) represents the expression value of each gene. The TCGA GC cohort was divided into low- and high-risk groups based on the median risk score, and then Kaplan–Meier analysis (R package: survival and survminer) was performed to analyze and compare the overall survival (OS) between the two subgroups. The time-dependent receiver operating characteristic (ROC) curve was adopted (R package: timeROC) to analyze the predictive accuracy and risk scores. Considering the clinical features, we also created a prognostic nomogram (R package: rms and survival) to predict the 1-, 3-, and 5-year overall survival.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) is a computational method that determines whether a pre-defined set of genes shows statistically significant differences between two biological states and is usually performed to estimate the changes in pathways and biological process activities in gene expression dataset samples. To assess the differences in biological processes between the high- and low-risk groups based on the gene expression profiling datasets, we downloaded the reference gene sets (c2. cp.v7.2. symbols.gmt and c5. all.v7.2. symbols.gmt) from the MSigDB and used the GSEA method included in the clusterProfiler R package (Yu et al., 2012) to perform enrichment analysis of the datasets and visualize them. An adjusted p-value of <0.05 and a false discovery rate (FDR, q-value) < 0.25 were considered to suggest statistically significant differences.

Estimation of immune cell-type fractions

The CIBERSORT (Chen et al., 2018) computational method was adopted to quantitatively convert the transcriptomic data of tumor tissues into the absolute abundance of immune cells and stromal cells, to evaluate the changes in the proportion of 22 human immune cell subpopulations. For each tumor sample, the sum of all the estimated immune cell-type fractions equaled 1.

Results

Analysis of 25 cell clusters from biopsy specimens of the human gastric antral mucosa with scRNA-seq data revealed high levels of cellular heterogeneity

After filtration based on the quality control criteria and normalization of the scRNA-seq data, 48,566 cells were obtained (Supplementary Figure S1A). We selected 2,000 highly variable genes for downstream analysis and tagged the top ten (Figure 1B), such as IGLL5, LIPF, TPSAB1, and APOA1. PCA was performed to identify available dimensions and relevant genes, and 20 principal components (PC) were selected for subsequent analysis (Figure 1C). The tSNE algorithm was successfully applied to divide human cells of the gastric antral mucosa into 25 independent clusters (Figure 1D), which were identified by the marker genes for each cell type (Supplementary Figures S1B,C). Table 1 and Supplementary Figure S1D showed the numbers and percentages of each cell type identified with the 25 clusters. Specifically, gastric epithelial, metaplastic stem-like, and enterocytes accounted for 55.53%, 15.59%, and 11%, of the total cell count, respectively. We then identified DEGs between cell types (Supplementary Figure S1E) and used the top two genes with the most significant differential expressions to draw dot plots (Figure 1E). We also calculated the number of each cell type and the proportion of cells per cell type in the samples of the different disease groups (Figures 1F,G). The aforementioned results portray the diverse landscape of the microenvironment between tumor and non-tumor samples.

FIGURE 1

Intercellular communication displayed locoregional immunomodulation in the carcinogenic process

We detected a total of 11 signaling pathways in the 11 cell types annotated in the single-cell data using CellChat, including MIF, MK, PTN, PARs, and GALECTIN. Heat maps were generated to illustrate the contribution of each pathway to the incoming or outgoing signals among the cell types (Figure 2A, Supplementary Figure S2). We selected the MIF signaling pathway as it contributed the most amongst the cell types, used a circle plot to illustrate the intensity of cell–cell interactions (Figure 2B) and performed a network centrality analysis (Figure 2C) for this pathway. Clearly, the macrophage is the most dominant recipient; B cell, T cell, and metaplastic stem-like cell perform key roles in the network. Moreover, we demonstrated the expression patterns of all ligand-receptor pairs included in the MIF signaling pathway (Figure 2D,E; Supplementary Figure S3) in different cell types: CD74 is mainly expressed in macrophages, endothelial cells, and B cells; CD44 is predominantly expressed in mast cells; and CXCR4 is primarily expressed in T cells.

FIGURE 2

Identification and enrichment analyses of differentially expressed genes between the early gastric cancer and IM groups

By analyzing the differences between the EGC and IM groups, we identified 64 DEGs in total (p-value < 0.05 and |log2FoldChange| > 1). Heat maps (Supplementary Figure S4) were drawn to visualize the DEGs, and a correlation analysis of these genes was performed (Supplementary Figure S5). As shown in Figures 3A,B, GKN1 and TFF2 are highly correlated in gastric epithelial cells; IGJ and IGLL5 are highly correlated in B cells; and OLFM4, REG1A, and TSPAN8 are highly correlated in gastric epithelial cells. GO and KEGG pathway enrichment analyses of the DEGs (adj. p value < 0.05) were visualized using dot plots (Figures 3C,D). These results prove that cell communication pathways (cadherin binding, adj. p = 2.5 × 10⁻⁴; ficolin−1−rich granule, adj. p= 8.6 × 10⁻⁶; regulation of cell−cell adhesion, adj. p = 7.0 × 10⁻⁵; and estrogen signaling pathways) and inflammatory response pathways (MHC protein complex, adj. p= 1.6 × 10⁻⁹; endocytic vesicle, adj. p = 2.6 × 10⁻⁵; antigen processing and presentation, and Th17 cell differentiation) were enriched.

FIGURE 3

Incidence of somatic mutations and copy number variations in 64 differentially expressed genes

We next examined the mutational landscape of 64 DEGs in TCGA patients with STAD. We found that the missense mutations have the highest mutation rate and the single-nucleotide polymorphism (SNP) is the most occurred variant type (Figure 4A). Among these SNPs, the highest proportion of mutations occurred in the C to T transition mutations. FCGBP and MUC6 are dominant in the top 10 mutated genes, and TMPRSS15, VIM, and APOA1 are also included. Figure 4B showed that the HSPH1 amplification (37%) and MUC6 (29%) and APOA1 (22%) frameshift deletions are the most common alterations observed. The mutagenesis in these TME-related DEGs is illustrated as the pro-carcinogenic potential.

FIGURE 4

Construction of a TME-related prognostic model based on prognostic genes

To evaluate the correlation between the prognosis for patients with GC and the expression of the 64 DEGs, we used univariate Cox regression analysis. The results suggested that six of the genes (TMPRSS15, VIM, LGALS1, APOA1, RNASE1, and TSC22D3) are significantly correlated with the disease prognosis (p < 0.05) (Figure 5A). In addition, the LASSO-Cox regression algorithm was used to establish a prognostic model, and the results suggested that four genes (TMPRSS15, VIM, APOA1, and RNASE1) are highly correlated with disease prognosis (Figure 5B). Risk score = (TMPRSS15^∗ 0.076) + (VIM^∗ 0.225) + (APOA1^∗ 0.066) + (RNASE1^∗ 0.135). The Kaplan–Meier plotter is an online tool to find survival biomarkers in GC based on the meta-analysis of GEO, TCGA, and European Genome-phenome Archive databases (Szász et al., 2016). Using the Kaplan–Meier plotter that has 875 enrolled patients with GC (Supplementary Figure S6), we further validated that these four genes have diagnostic value in association with adverse outcomes. We also verified their protein expression in the Human Protein Atlas database (Supplementary Figure S7), which is corroborated by the transcript levels in Supplementary Figure S4.

FIGURE 5

TCGA cohort was then divided into low- and high-risk groups based on the median risk score. The Kaplan–Meier curve suggested a lower overall survival rate for patients with GC with high-risk scores than those with low-risk scores (Figure 5C). The area under the curve was more significant than 0.6 for the 1-, 3-, and 5-year survival ROC curves, which indicated that risk scores accurately predicted patient survival rates (Figure 5D). The dataset’s risk score distribution and survival status are shown in Figures 4E,F, respectively. An increase in risk scores was accompanied by an increased risk of patient mortality and a shorter survival time (Figures 5E,F). Given the impact of clinicopathological features (including, age, gender, and TNM staging information) on prognosis, a prognostic nomogram was constructed to predict the survival of patients with GC in TCGA dataset (Figure 6A). The C-index of our nomogram model was 0.683, and the calibration plot (Figure 6B) showed a good fit with actual survival outcomes.

FIGURE 6

GSEA

GSEA was performed on all genes to analyze inter-group differences between the low- and high-risk groups. Figure 7 shows the nine most important functions or pathways based on normalized enrichment scores, such as GOBP: B cell receptor signaling pathway, GOCC: immunoglobulin complex circulating, and GOMF: antigen binding. As expected, increased risk scores are substantially associated with high-level immune response.

FIGURE 7

Immune cell infiltration

Based on the previous results, we found that the survival status of the high-risk group was significantly worse than that of the low-risk group and hence speculated that there might exist differences in immune cell infiltration between the low- and the high-risk groups. The CIBERSORT algorithm was used to calculate immune cell infiltration in GC. The bar graph of immune cell infiltration and the boxplot comparing immune cells between the low- and high-risk groups are shown in Figures 8A,B, respectively. We found differences between the two groups in M0 and M2 macrophages, resting mast cells, monocytes, resting NK cells, and CD8 and follicular helper T cells. Interestingly, the expression levels of M0 macrophages, resting NK cells, and follicular helper T cells increased in the low-risk group, which expressed a relatively mild immune response. In contrast, the expression levels of tumor-associated immune cells increased in the high-risk group such as M2 macrophages, resting mast cells, monocytes, and CD8 T cells. Considering that the RNA-seq data of TCGA referred to the expression levels of tissue blocks instead of blood samples, these tumor-immune infiltrates often indicate higher malignancy and worse prognosis rather than protective immunity.

FIGURE 8

Figure 9 shows the correlations between the expression of prognostic genes (TMPRSS15, VIM, APOA1, and RNASE1) and immune cell infiltration. APOA1 expression was positively correlated with the monocyte infiltration level; RNASE1 expression was positively correlated with the infiltration levels of M2 macrophages and CD8 T cells, and negatively correlated with the infiltration levels of M0 macrophages, resting NK cells, and follicular helper T cells; TMPRSS15 expression was negatively correlated with the M0 macrophages infiltration level; and VIM expression was positively correlated with the infiltration levels of M2 macrophages and resting mast cells, and negatively correlated with the follicular helper T cells infiltration level. This evidence suggested that the prognostic genes for GC are significantly associated with tumor immune infiltration.

FIGURE 9

Discussion

Among the numerous advances in the treatment of GC over the last decades, anti-PD-1/PD-L1 immunotherapy has raised particular concerns (Joshi and Badgwell, 2021). TME in GC may influence the response of immunotherapy and affect prognosis (Moutafi and Rimm, 2021). As transcriptomics technologies continue to advance, scRNA-seq offers an excellent tool to research the immune microenvironment biomarkers with clinical translational potential.

The first step toward single-cell analysis is to cluster cell populations with characteristic genes, which is the basis for mapping cell interactions. In addition to the common methods used in this study (Andrews and Hemberg, 2018), a new ensemble random projection-based algorithm, SHARP, which can cluster 10 million cells in large-scale scRNA-seq data, has recently appeared (Wan et al., 2020). CNVs analysis is usually applied for the identification of malignant tumor cells (Patel et al., 2014), but it has unsatisfactory effectiveness in GC (Zhang M. et al., 2021). Our research avoided the controversial topic of separating benign and malignant epithelial cells and instead categorized the component of TME. We found significant differences in the percentage of infiltrating immune cells between early GC and non-cancerous samples and between the low- and high-risk groups in the TCGA. These findings are consistent with a previous comparison of samples before and after cisplatin chemotherapy (Kim et al., 2021). To better clarify the humoral and cellular immune responses in GC, Sathe et al. (2020) concurrently sequenced matched peripheral blood mononuclear cells and revealed the immune remodeling of NK cells, dendritic cells, cytotoxic T cells, and plasma cells. Considering the loss of spatial dimension in cell isolation from tissues, the conjoint analysis of scRNA-seq and spatial transcriptomics in pancreatic cancer (Moncada et al., 2020) and squamous cell carcinoma (Ji et al., 2020) has been proposed, giving us a new perspective on tumor-infiltrating immune cells.

Some differentially expressed non-marker genes confirmed clustering results from another perspective. Overexpression of GKN1 in epithelial cells (Yoon et al., 2013), IGLL5 in B cells (Zheng et al., 2021), TPSAB1 in mast cells (Konnikova et al., 2021), and RGS5 in smooth muscle cells (Silini et al., 2012) are consistent with previous studies, while CXCL14 in fibroblasts (Wu et al., 2020) suggests that the conclusion from bulk RNA-seq may not be precise or versatile enough. Together with marker genes, these genes constitute 64 DEGs between IM and EGC, which reflect transcriptional heterogeneity among different cell types. In addition, 9 of 64 DEGs are members of the heat-shock protein family involved in cellular stress responses and pro-tumor inflammation as molecular chaperones. For example, Hsp72 has been previously reported to promote the oxaliplatin resistance of GC cells by inhibiting SDF-2 degradation (Takahashi et al., 2016); Hsp90ab1 is known to facilitate the epithelial–mesenchymal transition in GC by preventing LRP5 ubiquitination (Wang et al., 2019). Both genes are highly expressed by enterocytes and can be perceived as indicators of carcinogenicity for extensive intestinal metaplasia.

Macrophage migration inhibitory factor (MIF) is a multipotent cytokine involved in both inflammatory processes and anti-tumor immune response, having the properties of an enzyme, chemokine, and hormone simultaneously (Sumaiya et al., 2021). As a proinflammatory mediator secreted by numerous immune cells, MIF promotes inflammation and autoimmune diseases mainly by binding with the receptor CD74 and co-receptor CD44, CXCR4, or CXCR2 (Liehn et al., 2013; Nakamura et al., 2021; Wallace et al., 2021). Previous studies have shown that the carcinogenic role MIF plays is related to the activation of p53, Ras/MAPK, and Akt pathways in the inflammation–cancer axis (Mittelbronn et al., 2011; Wang et al., 2017; Liu et al., 2018). This study observed the prominent MIF signaling pathway network in GC at the single-cell level. Comparisons and enrichment analysis for DEGs in IM and EGC supported the enriched pathway associated with cell communication and inflammatory response. Mechanistically, MIF-pathway overactivation is likely a manifestation of either the macrophage aggregation caused by Helicobacter pylori infection or accumulations of tumor-associated macrophages, which further predisposes epithelial cells to malignant transformation. Under various environmental stimuli, macrophages developing from differentiated monocytes are classified into classically activated (M1) or alternatively activated (M2) macrophages. Generally, IL-12 synthesized and secreted by M1 can induce the proliferation and differentiation of naïve T cells into Th1 cells and enhance the NK cell-mediated antitumor effect; IL-10 secreted by M2 has the opposite effect (Gambardella et al., 2020). Our prognostic model shows that the expression of pro-tumorigenic M2 in the high-risk group is increased, and prognostic genes are also positively correlated with M2 infiltration. A reasonable explanation is that MIF facilitates the M2 polarization of macrophages in GC (Huang et al., 2019; Chen et al., 2021). Recently scholars (Zhou et al., 2022) found that carfilzomib enables M2 to express M1 cytokines, showing great immunotherapeutic potential for solid tumors. For future applications, experimental validation and subgroup analysis of macrophages are required.

Physician-scientists have long grappled with the quantitative assessment of immune TME (Zhang Z. et al., 2021). In advanced GC, Zeng et al. (2021) built a TME scoring system for the effects of checkpoint immunotherapy. By implementing machine learning, Cai et al. (2020) established a prognostic classification model based on the expression of VCAN, CLIP4, and MATN3 in a total of 1699 GC patients. In our study, four high-risk genes derived from scRNA-seq show the prognostication of poor outcomes. Vimentin (VIM), a type III intermediate filament protein, which characterizes the stromal component of TME in solid cancers, is often regarded as the epithelial–mesenchymal transition marker. It is closely related to GC cell invasion and metastasis, and is the most commonly used marker to detect the acquisition of these mesenchymal traits (Zeng et al., 2018; Zhu et al., 2019). Apolipoprotein A1 (APOA1), a component of high-density lipoprotein, exerts a favorable effect on the prevention of many cardiovascular diseases. However, its elevated level in the urine of bladder cancer patients (Peng et al., 2019) indicated shorter survival. For patients with colorectal cancer, increased APOA1 expression in the blood is concomitant with CD3⁺ T cells aggregation in the core of the tumor as well as the invasive margin (Guo et al., 2019). TMPRSS15, known as an enzyme gene, is translated as a serine protease in enterocytes and goblet cells—enteropeptidase. The link between TMPRSS15 and tumors has rarely been reported before, but other genes coding for the same family proteins were proven to be promising therapeutic targets. TMPRSS2 was found to favor the immune escape of COVID-19 (Shang et al., 2020), and TMPRSS4 was upregulated in malignancies of the stomach (Tazawa et al., 2022), liver (Dong et al., 2020), and prostate (Lee et al., 2021). In the present study, we speculate that TMPRSS15 and APOA1 used for enterocyte identification might be a cellular trait of intestinal type-GC. Therefore, the high expression of these two genes heralded extensive metaplasia and worse outcome. Similarly, despite the little evidence for a direct relationship between RNase1 release in endothelial cells (Bedenbender and Schmeck, 2020) and cancer progression, RNASE1 shows a good prediction performance based on raw data in TCGA (Li et al., 2020; Yang et al., 2022).

In conclusion, we analyzed the previously published data of scRNA-seq as well as TCGA, explored the inflammation–immunity–cancer axis, and developed a tumor microenvironment-associated risk score in GC. They well reflect the immune infiltration level and help construct a prognostic nomogram model that can assess overall long-term outcomes. Nevertheless, further in-depth studies are needed to confirm our results and broaden our view of TME in GC.

References

• 1

  AndrewsT. S.HembergM. (2018). Identifying cell populations with scRNASeq. Mol. Asp. Med.59, 114–122. 10.1016/j.mam.2017.07.002

  • CrossRef
  • Google Scholar

• 2

  BedenbenderK.SchmeckB. T. (2020). Endothelial ribonuclease 1 in cardiovascular and systemic inflammation. Front. Cell Dev. Biol.8, 576491. 10.3389/fcell.2020.576491

  • CrossRef
  • Google Scholar

• 3

  BijlsmaM.SadanandamA.TanP.VermeulenL. (2017). Molecular subtypes in cancers of the gastrointestinal tract. Nat. Rev. Gastroenterol. Hepatol.14 (6), 333–342. 10.1038/nrgastro.2017.33

  • CrossRef
  • Google Scholar

• 4

  CaiW. Y.DongZ. N.FuX. T.LinL. Y.WangL.YeG. D.et al (2020). Identification of a tumor microenvironment-relevant gene set-based prognostic signature and related therapy targets in gastric cancer. Theranostics10 (19), 8633–8647. 10.7150/thno.47938

  • CrossRef
  • Google Scholar

• 5

  ChenB.KhodadoustM. S.LiuC. L.NewmanA. M.AlizadehA. A. (2018). Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol. Biol.1711, 243–259. 10.1007/978-1-4939-7493-1_12

  • CrossRef
  • Google Scholar

• 6

  ChenJ.WangH.JiaL.HeJ.LiY.LiuH.et al (2021). Bufalin targets the SRC-3/MIF pathway in chemoresistant cells to regulate M2 macrophage polarization in colorectal cancer. Cancer Lett.513, 63–74. 10.1016/j.canlet.2021.05.008

  • CrossRef
  • Google Scholar

• 7

  DongZ. R.SunD.YangY. F.ZhouW.WuR.WangX. W.et al (2020). TMPRSS4 drives angiogenesis in hepatocellular carcinoma by promoting HB-egf expression and proteolytic cleavage. Hepatology72 (3), 923–939. 10.1002/hep.31076

  • CrossRef
  • Google Scholar

• 8

  GambardellaV.CastilloJ.TarazonaN.Gimeno-ValienteF.Martinez-CiarpagliniC.Cabeza-SeguraM.et al (2020). The role of tumor-associated macrophages in gastric cancer development and their potential as a therapeutic target. Cancer Treat. Rev.86, 102015. 10.1016/j.ctrv.2020.102015

  • CrossRef
  • Google Scholar

• 9

  GuoG.WangY.ZhouY.QuanQ.ZhangY.WangH.et al (2019). Immune cell concentrations among the primary tumor microenvironment in colorectal cancer patients predicted by clinicopathologic characteristics and blood indexes. J. Immunother. Cancer7 (1), 179. 10.1186/s40425-019-0656-3

  • CrossRef
  • Google Scholar

• 10

  HanahanD. (2022). Hallmarks of cancer: New dimensions. Cancer Discov.12 (1), 31–46. 10.1158/2159-8290.CD-21-1059

  • CrossRef
  • Google Scholar

• 11

  HuangW. C.KuoK. T.WangC. H.YehC. T.WangY. (2019). Cisplatin resistant lung cancer cells promoted M2 polarization of tumor-associated macrophages via the Src/CD155/MIF functional pathway. J. Exp. Clin. Cancer Res.38 (1), 180. 10.1186/s13046-019-1166-3

  • CrossRef
  • Google Scholar

• 12

  JeongH. Y.HamI. H.LeeS. H.RyuD.SonS. Y.HanS. U.et al (2021). Spatially distinct reprogramming of the tumor microenvironment based on tumor invasion in diffuse-type gastric cancers. Clin. Cancer Res.27 (23), 6529–6542. 10.1158/1078-0432.CCR-21-0792

  • CrossRef
  • Google Scholar

• 13

  JiA. L.RubinA. J.ThraneK.JiangS.ReynoldsD. L.MeyersR. M.et al (2020). Multimodal analysis of composition and spatial architecture in human squamous cell carcinoma. Cell182 (2), 1661–1662. e422. 10.1016/j.cell.2020.08.043

  • CrossRef
  • Google Scholar

• 14

  JinS.Guerrero-JuarezC. F.ZhangL.ChangI.RamosR.KuanC. H.et al (2021). Inference and analysis of cell-cell communication using CellChat. Nat. Commun.12 (1), 1088. 10.1038/s41467-021-21246-9

  • CrossRef
  • Google Scholar

• 15

  JoshiS. S.BadgwellB. D. (2021). Current treatment and recent progress in gastric cancer. Ca. Cancer J. Clin.71 (3), 264–279. 10.3322/caac.21657

  • CrossRef
  • Google Scholar

• 16

  KimR.AnM.LeeH.MehtaA.HeoY. J.KimK. M.et al (2021). Early tumor–immune microenvironmental remodeling and response to first-line fluoropyrimidine and platinum chemotherapy in advanced gastric cancer. Cancer Discov.12, 984–1001. 10.1158/2159-8290.CD-21-0888

  • CrossRef
  • Google Scholar

• 17

  KonnikovaL.RobinsonT. O.OwingsA. H.ShirleyJ. F.DavisE.TangY.et al (2021). Small intestinal immunopathology and GI-associated antibody formation in hereditary alpha-tryptasemia. J. Allergy Clin. Immunol.148 (3), 813–821.e7. e817. 10.1016/j.jaci.2021.04.004

  • CrossRef
  • Google Scholar

• 18

  KumarV.RamnarayananK.SundarR.PadmanabhanN.SrivastavaS.KoiwaM.et al (2021). Single-cell Atlas of lineage states, tumor microenvironment, and subtype-specific expression programs in gastric cancer. Cancer Discov.12, 670–691. 10.1158/2159-8290.CD-21-0683

  • CrossRef
  • Google Scholar

• 19

  LeeY.YoonJ.KoD.YuM.LeeS.KimS.et al (2021). TMPRSS4 promotes cancer stem-like properties in prostate cancer cells through upregulation of SOX2 by SLUG and TWIST1. J. Exp. Clin. Cancer Res.40 (1), 372. 10.1186/s13046-021-02147-7

  • CrossRef
  • Google Scholar

• 20

  LiS.JiangL.TangJ.GaoN.GuoF. (2020). Kernel fusion method for detecting cancer subtypes via selecting relevant expression data. Front. Genet.11, 979. 10.3389/fgene.2020.00979

  • CrossRef
  • Google Scholar

• 21

  LiehnE.KanzlerI.KonschallaS.KrohA.SimsekyilmazS.SönmezT.et al (2013). Compartmentalized protective and detrimental effects of endogenous macrophage migration-inhibitory factor mediated by CXCR2 in a mouse model of myocardial ischemia/reperfusion. Arterioscler. Thromb. Vasc. Biol.33 (9), 2180–2186. 10.1161/atvbaha.113.301633

  • CrossRef
  • Google Scholar

• 22

  LiuR.SunD.JiaoY.WangP.ZhangJ.WangM.et al (2018). Macrophage migration inhibitory factor promotes tumor aggressiveness of esophageal squamous cell carcinoma via activation of Akt and inactivation of GSK3β. Cancer Lett.412, 289–296. 10.1016/j.canlet.2017.10.018

  • CrossRef
  • Google Scholar

• 23

  MayakondaA.LinD. C.AssenovY.PlassC.KoefflerH. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res.28 (11), 1747–1756. 10.1101/gr.239244.118

  • CrossRef
  • Google Scholar

• 24

  MittelbronnM.PlattenM.ZeinerP.DombrowskiY.FrankB.ZachskornC.et al (2011). Macrophage migration inhibitory factor (MIF) expression in human malignant gliomas contributes to immune escape and tumour progression. Acta Neuropathol.122 (3), 353–365. 10.1007/s00401-011-0858-3

  • CrossRef
  • Google Scholar

• 25

  MoncadaR.BarkleyD.WagnerF.ChiodinM.DevlinJ. C.BaronM.et al (2020). Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat. Biotechnol.38 (3), 333–342. 10.1038/s41587-019-0392-8

  • CrossRef
  • Google Scholar

• 26

  MoutafiM.RimmD. (2021). Putting the microenvironment into the immunotherapy companion diagnostic. Clin. Cancer Res.27 (14), 3812–3814. 10.1158/1078-0432.ccr-21-1238

  • CrossRef
  • Google Scholar

• 27

  MullerS.DiazA. (2017). Single-cell mRNA sequencing in cancer research: Integrating the genomic fingerprint. Front. Genet.8, 73. 10.3389/fgene.2017.00073

  • CrossRef
  • Google Scholar

• 28

  MundekkadD.ChoW. C. (2022). Nanoparticles in clinical translation for cancer therapy. Int. J. Mol. Sci.23 (3), 1685. 10.3390/ijms23031685

  • CrossRef
  • Google Scholar

• 29

  NakamuraA.ZengF.NakamuraS.ReidK.GraceyE.LimM.et al (2021). Macrophage migration inhibitory factor drives pathology in a mouse model of spondyloarthritis and is associated with human disease. Sci. Transl. Med.13 (616), eabg1210. 10.1126/scitranslmed.abg1210

  • CrossRef
  • Google Scholar

• 30

  PatelA. P.TiroshI.TrombettaJ. J.ShalekA. K.GillespieS. M.WakimotoH.et al (2014). Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science344 (6190), 1396–1401. 10.1126/science.1254257

  • CrossRef
  • Google Scholar

• 31

  PengC.HuaM. Y.LiN. S.HsuY. P.ChenY. T.ChuangC. K.et al (2019). A colorimetric immunosensor based on self-linkable dual-nanozyme for ultrasensitive bladder cancer diagnosis and prognosis monitoring. Biosens. Bioelectron.126, 581–589. 10.1016/j.bios.2018.11.022

  • CrossRef
  • Google Scholar

• 32

  PicelliS.FaridaniO. R.BjorklundA. K.WinbergG.SagasserS.SandbergR.et al (2014). Full-length RNA-seq from single cells using Smart-seq2. Nat. Protoc.9 (1), 171–181. 10.1038/nprot.2014.006

  • CrossRef
  • Google Scholar

• 33

  RamskoldD.LuoS.WangY. C.LiR.DengQ.FaridaniO. R.et al (2012). Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells. Nat. Biotechnol.30 (8), 777–782. 10.1038/nbt.2282

  • CrossRef
  • Google Scholar

• 34

  SatheA.GrimesS. M.LauB. T.ChenJ.SuarezC.HuangR. J.et al (2020). Single-cell genomic characterization reveals the cellular reprogramming of the gastric tumor microenvironment. Clin. Cancer Res.26 (11), 2640–2653. 10.1158/1078-0432.CCR-19-3231

  • CrossRef
  • Google Scholar

• 35

  SatijaR.FarrellJ. A.GennertD.SchierA. F.RegevA. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol.33 (5), 495–502. 10.1038/nbt.3192

  • CrossRef
  • Google Scholar

• 36

  ShangJ.WanY.LuoC.YeG.GengQ.AuerbachA.et al (2020). Cell entry mechanisms of SARS-CoV-2. Proc. Natl. Acad. Sci. U. S. A.117 (21), 11727–11734. 10.1073/pnas.2003138117

  • CrossRef
  • Google Scholar

• 37

  SiliniA.GhilardiC.FiginiS.SangalliF.FruscioR.DahseR.et al (2012). Regulator of G-protein signaling 5 (RGS5) protein: A novel marker of cancer vasculature elicited and sustained by the tumor's proangiogenic microenvironment. Cell. Mol. Life Sci.69 (7), 1167–1178. 10.1007/s00018-011-0862-8

  • CrossRef
  • Google Scholar

• 38

  SmythE. C.NilssonM.GrabschH. I.van GriekenN. C. T.LordickF. (2020). Gastric cancer. Lancet396 (10251), 635–648. 10.1016/s0140-6736(20)31288-5

  • CrossRef
  • Google Scholar

• 39

  SumaiyaK.LangfordD.NatarajaseenivasanK.ShanmughapriyaS. (2021). Macrophage migration inhibitory factor (MIF): A multifaceted cytokine regulated by genetic and physiological strategies. Pharmacol. Ther.233, 108024. 10.1016/j.pharmthera.2021.108024

  • CrossRef
  • Google Scholar

• 40

  SzászA.LánczkyA.NagyÁ.FörsterS.HarkK.GreenJ.et al (2016). Cross-validation of survival associated biomarkers in gastric cancer using transcriptomic data of 1, 065 patients. Oncotarget7 (31), 49322–49333. 10.18632/oncotarget.10337

  • CrossRef
  • Google Scholar

• 41

  TakahashiK.TanakaM.YashiroM.MatsumotoM.OhtsukaA.NakayamaK. I.et al (2016). Protection of stromal cell-derived factor 2 by heat shock protein 72 prevents oxaliplatin-induced cell death in oxaliplatin-resistant human gastric cancer cells. Cancer Lett.378 (1), 8–15. 10.1016/j.canlet.2016.05.002

  • CrossRef
  • Google Scholar

• 42

  TangF.BarbacioruC.WangY.NordmanE.LeeC.XuN.et al (2009). mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods6 (5), 377–382. 10.1038/nmeth.1315

  • CrossRef
  • Google Scholar

• 43

  TazawaH.SuzukiT.SaitoA.IshikawaA.KomoT.SadaH.et al (2022). Utility of TMPRSS4 as a prognostic biomarker and potential therapeutic target in patients with gastric cancer. J. Gastrointest. Surg.26 (2), 305–313. 10.1007/s11605-021-05101-2

  • CrossRef
  • Google Scholar

• 44

  WallaceD.FiguerasF.WegenerW.GoldenbergD. (2021). Experience with milatuzumab, an anti-CD74 antibody against immunomodulatory macrophage migration inhibitory factor (MIF) receptor, for systemic lupus erythematosus (SLE). Ann. Rheum. Dis.80 (7), 954–955. 10.1136/annrheumdis-2020-219803

  • CrossRef
  • Google Scholar

• 45

  WanS.KimJ.WonK. J. (2020). SHARP: Hyperfast and accurate processing of single-cell RNA-seq data via ensemble random projection. Genome Res.30 (2), 205–213. 10.1101/gr.254557.119

  • CrossRef
  • Google Scholar

• 46

  WangB.ZhangY.QingT.XingK.LiJ.ZhenT.et al (2021). Comprehensive analysis of metastatic gastric cancer tumour cells using single-cell RNA-seq. Sci. Rep.11 (1), 1141. 10.1038/s41598-020-80881-2

  • CrossRef
  • Google Scholar

• 47

  WangC.ZhouX.LiW.LiM.TuT.BaX.et al (2017). Macrophage migration inhibitory factor promotes osteosarcoma growth and lung metastasis through activating the RAS/MAPK pathway. Cancer Lett.403, 271–279. 10.1016/j.canlet.2017.06.011

  • CrossRef
  • Google Scholar

• 48

  WangH.DengG.AiM.XuZ.MouT.YuJ.et al (2019). Hsp90ab1 stabilizes LRP5 to promote epithelial-mesenchymal transition via activating of AKT and Wnt/β-catenin signaling pathways in gastric cancer progression.Oncogene38 (9), 1489–1507. 10.1038/s41388-018-0532-5

  • CrossRef
  • Google Scholar

• 49

  WangR.DangM.HaradaK.HanG.WangF.Pool PizziM.et al (2021). Single-cell dissection of intratumoral heterogeneity and lineage diversity in metastatic gastric adenocarcinoma. Nat. Med.27 (1), 141–151. 10.1038/s41591-020-1125-8

  • CrossRef
  • Google Scholar

• 50

  WangR.ZhengX.WangJ.WanS.SongF.WongM. H.et al (2022). Improving bulk RNA-seq classification by transferring gene signature from single cells in acute myeloid leukemia. Brief. Bioinform.23 (2), bbac002. 10.1093/bib/bbac002

  • CrossRef
  • Google Scholar

• 51

  WuD.LiZ.ZhaoS.YangB.LiuZ. (2020). Downregulated microRNA-150 upregulates IRX1 to depress proliferation, migration, and invasion, but boost apoptosis of gastric cancer cells. IUBMB Life72 (3), 476–491. 10.1002/iub.2214

  • CrossRef
  • Google Scholar

• 52

  YangH.TianW.ZhouB. (2022). Sarcopenia and a 5-mRNA risk module as a combined factor to predict prognosis for patients with stomach adenocarcinoma. Genomics114 (1), 361–377. 10.1016/j.ygeno.2021.12.011

  • CrossRef
  • Google Scholar

• 53

  YoonJ. H.ChoiY. J.ChoiW. S.AshktorabH.SmootD. T.NamS. W.et al (2013). GKN1-miR-185-DNMT1 axis suppresses gastric carcinogenesis through regulation of epigenetic alteration and cell cycle. Clin. Cancer Res.19 (17), 4599–4610. 10.1158/1078-0432.CCR-12-3675

  • CrossRef
  • Google Scholar

• 54

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 55

  ZengD.WuJ.LuoH.LiY.XiaoJ.PengJ.et al (2021). Tumor microenvironment evaluation promotes precise checkpoint immunotherapy of advanced gastric cancer. J. Immunother. Cancer9 (8), e002467. 10.1136/jitc-2021-002467

  • CrossRef
  • Google Scholar

• 56

  ZengS.XieX.XiaoY. F.TangB.HuC. J.WangS. M.et al (2018). Long noncoding RNA LINC00675 enhances phosphorylation of vimentin on Ser83 to suppress gastric cancer progression. Cancer Lett.412, 179–187. 10.1016/j.canlet.2017.10.026

  • CrossRef
  • Google Scholar

• 57

  ZengY.JinR. (2021). Molecular pathogenesis, targeted therapies, and future perspectives for gastric cancer. Semin. Cancer Biol2021, 1044–579X. 10.1016/j.semcancer.2021.12.004

  • CrossRef
  • Google Scholar

• 58

  ZhangM.HuS.MinM.NiY.LuZ.SunX.et al (2021). Dissecting transcriptional heterogeneity in primary gastric adenocarcinoma by single cell RNA sequencing. Gut70 (3), 464–475. 10.1136/gutjnl-2019-320368

  • CrossRef
  • Google Scholar

• 59

  ZhangP.YangM.ZhangY.XiaoS.LaiX.TanA.et al (2019). Dissecting the single-cell transcriptome network underlying gastric premalignant lesions and early gastric cancer. Cell Rep.27 (6), 1934–1947. 10.1016/j.celrep.2019.04.052

  • CrossRef
  • Google Scholar

• 60

  ZhangZ.BaoS.YanC.HouP.ZhouM.SunJ.et al (2021). Computational principles and practice for decoding immune contexture in the tumor microenvironment. Brief. Bioinform.22 (3), bbaa075. 10.1093/bib/bbaa075

  • CrossRef
  • Google Scholar

• 61

  ZhaoL.LeeV. H. F.NgM. K.YanH.BijlsmaM. F. (2019). Molecular subtyping of cancer: Current status and moving toward clinical applications. Brief. Bioinform.20 (2), 572–584. 10.1093/bib/bby026

  • CrossRef
  • Google Scholar

• 62

  ZhengB.XieF.ChengF.WangJ.YaoZ.HeW.et al (2021). Integrative analysis of immune-related genes in the tumor microenvironment of renal clear cell carcinoma and renal papillary cell carcinoma. Front. Mol. Biosci.8, 760031. 10.3389/fmolb.2021.760031

  • CrossRef
  • Google Scholar

• 63

  ZhouQ.LiangJ.YangT.LiuJ.LiB.LiY.et al (2022). Carfilzomib modulates tumor microenvironment to potentiate immune checkpoint therapy for cancer. EMBO Mol. Med.14 (1), e14502. 10.15252/emmm.202114502

  • CrossRef
  • Google Scholar

• 64

  ZhuZ.RongZ.LuoZ.YuZ.ZhangJ.QiuZ.et al (2019). Circular RNA circNHSL1 promotes gastric cancer progression through the miR-1306-3p/SIX1/vimentin axis. Mol. Cancer18 (1), 126. 10.1186/s12943-019-1054-7

  • CrossRef
  • Google Scholar