Identification of 10 Gene Signatures as Prognostic Biomarkers for Peritoneal Metastatic Gastric Cancer

Abstract

Background. Peritoneal metastatic gastric cancer (PMGC) is very common, and usually, the prognosis is poor. There is currently an absence of accurate methods for the early diagnosis and prediction of peritoneal metastasis (PM). This highlights the need to develop strategies to identify the risk of PMGC.

Methods. We performed a comprehensive discovery of biomarkers to predict PM by analyzing profiling datasets from GSE62254. The prognostic PM-related genes were obtained using the univariate Cox regression analysis, followed by a least absolute shrinkage and selection operator regression (LASSO) to establish a risk score model. The gene set enrichment analysis (GSEA) was used to determine the pathway enrichment in both the high- and low-risk groups. The 1-, 3-, and 5-year overall survival (OS) rates and area under the receiver operating characteristic curve (ROC) were used to compare the predictive accuracy-based risk stratification. In addition, an unsupervised clustering algorithm was applied to divide patients into subgroups according to the PM-related genes.

Results. We identified 10 genes (MMP12, TAC1, TSPYL5, PPP1R14A, TMSB15B, NPY1R, PCDH9, EPM2AIP1, TIG7, and DYNC1I1) for PMGC diagnosis. The OS rates between the high- and low-risk groups at 1-, 3-, and 5-years were significantly different in the training and validation sets. The AUCs at 1-, 3-, and 5-years in the training set were 0.71, 0.74, and 0.73, respectively. In the validation set, the AUCs at 1-, 3-, and 5-years were 0.68, 0.66, and 0.69, respectively. The 10 gene signatures were correlated with immune cell infiltration in both the high- and low-risk groups. In addition, based on the GSEA, several significant pathways were enriched in the high-risk PMGC group, such as the Wnt and transforming growth factor beta (TGF-β) signaling pathway and leukocyte transendothelial migration pathway. Furthermore, unsupervised cluster analysis showed that the model could distinguish the level of risk among patients with PMGC.

Conclusions. Overall, 10 gene signatures were identified for PMGC risk prediction. These may be valuable in making clinical decisions to improve treatment outcomes in patients with PMGC.

Introduction

Gastric cancer (GC) has the sixth highest incidence rate and is the second leading cause of cancer-related deaths globally, with7.7% of total cancer-related deaths in 2020[1]. Approximately 70% of the cases are diagnosed in developing countries, half of which are in East Asia[2]. Although the comprehensive treatment for gastric cancer has improved, the 5-year overall survival (OS) rate for advanced gastric cancer is only around 31–69% due to its high heterogeneity[3]. Moreover, gastric cancer is commonly diagnosed during the end stages in those countries without national screening programs[4]. Among the different dissemination patterns, peritoneal metastasis (PM) is the most common, accounting for 53–66% of distant metastases in gastric cancer[5].

Poor prognosis is common among patients with PM. Even patients with positive peritoneal lavage cytology have a poor prognosis after radical gastrectomy[6, 7]. Furthermore, the response rate of initially diagnosed patients with PM is approximately 14%, while the median survival time is approximately 3–7 months, and the 5-year OS rate is 6%[8]. The presence of PM is also considered a contraindication to surgery, and instead, systemic adjuvant chemotherapy is recommended. Early diagnosis of PM is a challenge because of the lack of typical symptoms and difficulty in detecting small peritoneal deposits. The initial clinical manifestation of abdominal metastasis is ascites or intestinal obstruction[9]. Clinically, PM is diagnosed by laparoscopic exploration, abdominal lavage, and computed tomography[10, 11]; however, these methods have certain limitations. Laparoscopic exploration is invasive and entails additional risks, while peritoneal lavage cytological examination and CT have poor sensitivity when it comes to detecting intraperitoneal cancer cells. It has been reported that PM and non-PM subgroups in patients with GC have variations in their PM-related genes[12]. However, the substantial heterogeneity of GC depends on complicated pathway changes rather than changes in a single gene. Hence, it is crucial to distinguish subtype-specific multi-gene signatures for the prognostication and identification of molecular differences in peritoneal metastatic gastric cancer (PMGC). Lastly, the tumor microenvironment (TME)—a comprehensive system comprising tumor and stromal tissue, surrounded by tumor-related fibroblasts, inflammatory, immune cells and interstitial tissue, and various cytokines and chemokines, influences how a tumor grows and spreads. As such, analysis of immune cell infiltration in cancer tissue plays a crucial role in disease research and prediction of treatment prognosis[13].

In this study, we explored the gastric cancer genome-wide expression profile dataset GSE62254 via the public Gene Expression Omnibus (GEO) database to establish and to verify novel gene signatures for PMGC. Determining the accuracy of gene signatures in predicting the risk and OS rate among patients may help understand the potential prognostic value of PM-related genes in GC.

Materials & Methods

Data collection and differentially expressed gene (DEG) analysis

We analyzed the gastric cancer genome-wide expression profile data GSE62254[14] from the public GEO database[15] using the R package "GEOquery"[16]. GSE62254 contained 300 gastric cancer tissue cases, wherein 246 were non-PM and 54 were PM. We used the "limma" package[17] to explore the gene expression level differences between the PM and non-PM group to obtain DEGs. We set genes with logFC>0.1 at an adjusted p<0.05, as upregulated genes, and genes with logFC<0.1 at an adjusted p<0.05 as downregulated genes.

Gene set enrichment analysis (GSEA)

GSEA was utilized to assess the distribution tendency of genes. We used a predefined gene set in the gene list ranked by the phenotype correlation to determine its contribution to the phenotype (Subramanian et al., 2005). Specifically, the gene set C2.cp.kegg.v7.0.entrez.gmt was used in the GSEA in the data set GSE62254 using the R package "clusterProfiler" [18]with false discovery rate (FDR)<0.25. Statistical significance was set at a p-value of <0.05. We then evaluated the potential functions of the signature in dataset GSE62254 in comparing the high- and low-risk groups.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes(KEGG) analysis

Using the R package "clusterProfiler," we performed GO function annotation analysis and KEGG pathway enrichment analysis on DEGs related to PMGC to explore the potential functions of DEGs. Statistical significance was set at p<0.05.

Univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression model

Univariate regression was used to screen DEGs for prognosis-related PMGC genes. Three hundred patients were divided into training and validation sets with a 7:3 ratio. LASSO regression uses shrinkage to analyze the genes that meet the outcome of the Cox regression (p<0.0001) for further dimension reduction. Critical prognostic genes related to PM were screened to establish a LASSO regression risk model.

Single sample GSEA (ssGSEA) immune infiltration analysis

We used thegene set variation analysis package[19] to perform ssGSEA analysis of 300 cases of gastric cancer. This was performed to estimate the distribution and abundance of 27 immuno-related cells in the TME, to analyze the correlation between immune cells, and to compare the immune cell infiltration between the high- and low-risk groups.

Unsupervised clustering analysis

The unsupervised clustering algorithm was applied to divide patients into subgroups according to the essential peritoneal dissemination genes using the R package "ConsensusClusterPlus"[20]. Survival curves were drawn to compare patient survival. Statistical significance was set at p<0.05.

Statistical analyses

Statistical analyses were performed using R software version 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria)[21].The continuous variables of the two groups were compared by the Mann-Whitney U test. The differences and correlations determined with the Kruskal–Wallis test variables were analyzed using the Spearman algorithm. Statistical significance was set at p<0.05.

Results

Identification of DEGs

The workflow of this study is shown in Fig. 1. Overall, 300 patients from the GSE62254 dataset were divided into PM (n=54) and non-PM groups (n=246). The predictive analysis revealed significant differences in DEGs between the PM and non-PM groups. A total of 987 DEGs were identified, including 701 upregulated genes and 286 downregulated genes. Volcano maps and heat maps were generated for the DEGs (Fig. 2).

GSEA of potential pathways

Six pathways were enriched in the PMGC group: specifically, signaling pathways, such as the Wnt, TGF-β, calcium, and mitogen-activated protein kinase, along with the extracellular matrix-receptor interaction pathway, and gap junction (Fig. 3).

GO annotation and KEGG pathway enrichment

DEGs were enriched in the TGF-β and interleukin (IL)-17 signaling pathways, cell cycle, vascular smooth muscle contraction, and other biological pathways. The HSA04350: TGF-β and HSA04657: IL-17 signaling pathways were significantly enriched and are shown in Fig. 4.

Univariate Cox and LASSO regression

We used univariate Cox regression analysis for the 987 DEGs and obtained 116 prognostic PM-related genes (Table 1). LASSO regression screening resulted in 10 key genes and a LASSO–Cox regression model was established (Fig. 5). The expression value of each gene was calculated using the following formula:

Based on risk score, the training set was divided into two groups, high-risk group and low-risk group. The OS rate was significantly higher in the low-risk group than in the high-risk group in both the training (p=1.321×10^-8) and validation (p=0.0094) sets. The time-dependent AUCs under the receiver operating characteristic curve at 1-, 3-, and 5-years for the training set were 0.71, 0.74, and 0.73, respectively. In the validation set, the AUCs at 1-, 3-, and 5-years were 0.68, 0.66, and 0.69, respectively. The distribution of patients' risk scores and survival times are shown in the test and validation sets (Fig. 5).

GSEA of the high- and low-risk groups

GSEA results showed that tissue migration, negative regulation of leukocyte chemotaxis, negative regulation of leukocyte migration, endothelial cell migration, and macrophage cytokine production were all enriched in the high-risk group. In the low-risk group, RNA modification, somatic recombination of immunoglobulin gene segments, mitochondrial gene expression, and mitochondrial matrix were enriched. The results of the KEGG analysis showed that the Wnt andTGF-β signaling pathways andleukocyte transendothelial migration pathway were significantly enriched in the high-risk group. Meanwhile, the apoptosis and natural killer cell-mediated cytotoxicity pathways were significantly enriched in the low-risk group (Fig 6).

Immunoinfiltration analysis

The correlation heatmap showed correlations within 27 types of immune cells. Between the high- and low-risk groups, there were significant differences in the expression of 18 types of immune cells (Fig 7).

Unsupervised clustering analysis

Consensus clustering analysis was performed for 300 patients based on 10 prognostic PM genes (Fig. 8). Partition-assisted clusteringwas used to calculate the optimal clustering (K=4). Heatmaps were used to display the expression of 10 prognostic PM genes in four subcategories. The 10 molecules were significantly different among the four prognostic subgroups (Cluster 1, Cluster 2, Cluster 3, and Cluster 4). In addition, survival analysis showed significant differences in the OS rates of patients in the different subgroups (p=1.66 × 10^-11).

Discussion

In this study, we identified 116 prognostic PM-related genesand followed univariate and LASSO regression to establish a prognostic model of 10 genes (MMP12, TAC1, TSPYL5, PPP1R14A, TMSB15B, NPY1R, PCDH9, EPM2AIP1, TIG7, and DYNC1I1). These models showed good predictive values in the training and validation sets. We also compared the infiltrate patterns of immune cells between high- and low-risk groups of PMGC. Furthermore, based on the unsupervised cluster analysis, the model was able to effectively distinguish patients with different risk levels among those with PMGC.

Evidence accumulated in recent years has led to a better understanding of the various molecular subtypes of GC[22]. Microsatellitestable and epithelial‑mesenchymal transition GC subtypes have been found to have a higher likelihood of PM (64%) than the other three subtypes, indicating a high correlation between this subtype and PM[14]. In addition, in a multi-omic analysis of isolated corresponding tumor cell lines from malignant ascitic fluid samples from patients with GC, two molecular subtypes were identified to be possible clinical therapeutic targets[23] . The molecular complex method was also previously used to identify key genes and pathways for predicting PM[12]. Despite considerable efforts to uncover the genomic and epigenomic features of PMGC, its implications on the clinical disease management process for patients with GC, especially those with PMGC, remain limited. Data on molecular biomarkers that can predict peritoneal recurrence are still lacking. This is highly important for improving the disease management and treatment decisions for GC.

Recently, many studies have shown that immune cell infiltration in GC significantly impacts clinical outcomes[24-26]. Hence, we compared the different infiltrate immune cells related to GC between the high- and low-risk groups and noted 18 significantly different immune cells. Compared with the low-risk group, activated B cells, memory T cell (mTc), effector memory CD4 T cells, eosinophils, immature B cells, mast cells, memory B cells, monocytes, natural killer cells, plasmacytoid dendritic cells, T follicular helper cells, and type 1 T helper cells were significantly higher in the high-risk group. Comprehensive research has shown that plasmacytoid dendritic cells are correlated with poor prognosis in GC[27]. In renal cell carcinoma, a lower neutrophil-to-eosinophil ratio is predictive of betterclinical outcomes[28]. We also observed a high level of mTc in the high-risk group, and there is evidence that high mTc is correlated with tumor load[29]. Therefore, the abovementioned immune cell infiltration patterns may better explain why the high-risk group had a worse prognosis. Furthermore, we used GSEA to investigate the different biological processes between the two groups. The results revealed that six pathwayswere significantly enriched in the high-risk group. The TGF-β pathway has been shown to promote vascular endothelial growth factor production and to induce differentiation of cancer-associated fibroblasts and affect T-cells balance[30]. Some pathological and physiological processes, including cancer, are controlled by canonical and non-canonical Wnt/β-catenin pathways[31]. Thus, the identified enriched pathways are mainly linked to gastric peritoneal tumorigenesis, and further work to determine potential mechanisms will help pinpointmoreeffective therapeutic targets for GC.

We recognize that the current study has limitations. First, the data analysis was based on genomics rather than proteomics. Studies have shown that there are higher correlations between strains of quantitative RNA sequencing and proteomics data; however, these correlate poorly among themselves. This poor correlation may be extensive and may affect the accuracy of prediction models[32, 33]. Second, the prediction model was not validated using different independent datasets. Lastly, 10 genes would still be considered numerous for a prediction model and can affect its clinical application to some extent. Further research is still needed to include more validation studies along with utilizing machine learning to continue to improve the model.

Conclusions

In conclusion, our systematic and comprehensive biomarker discovery and validation work has enabled the identification of 10 genetic biomarkers that predict PMGC. This model can identify patients at a high risk for developing PMGC and may contribute to making appropriate clinical decisions and improving treatment outcomes in patients with GC.

Abbreviations

Declarations

Acknowledgements

None

Authors’ contributions

Guo Tiankang designed the study. Li Junliang and Zhang Lingfang searched and analysised the data. Li Junliang and Zhang Lingfang wrote the manuscript. All author read and approved the final manuscript.

Funding

This study was supported by Natural Science Foundation of Gansu Province (No. 20JR10RA372),Health Industry Research Project of Gansu Province (GSWSKY-2019-03),National scientific research project support program of Gansu Provincial Hospital (19SYPYB-10),Open fund project of Gansu Key Laboratory of Molecular Diagnosis and Precision Therapy of Surgical Oncology (2020GSZDSYS02), Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences (NLDTG2020020).

Availability of data and materials

The datasets analyzed were acquired from GEOdatabase(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62254).

Ethics approval and consent to participate

There was not applicable to this manuscript.

Consent for publication

Consent for publication was obtained from all participants.

Competing interests

The authors declare that they have no competing interests.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F: Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries.CA Cancer J Clin 2021, 71:209-249.
2. Rahman R, Asombang AW, Ibdah JA: Characteristics of gastric cancer in Asia.World J Gastroenterol 2014, 20:4483-4490.
3. Stahl KA, Olecki EJ, Dixon ME, Peng JS, Torres MB, Gusani NJ, Shen C: Gastric Cancer Treatments and Survival Trends in the United States.Curr Oncol 2020, 28:138-151.
4. Kim YG, Kong SH, Oh SY, Lee KG, Suh YS, Yang JY, Choi J, Kim SG, Kim JS, Kim WH, et al: Effects of screening on gastric cancer management: comparative analysis of the results in 2006 and in 2011.J Gastric Cancer 2014, 14:129-134.
5. Dong D, Tang L, Li ZY, Fang MJ, Gao JB, Shan XH, Ying XJ, Sun YS, Fu J, Wang XX, et al: Development and validation of an individualized nomogram to identify occult peritoneal metastasis in patients with advanced gastric cancer.Ann Oncol 2019, 30:431-438.
6. Higaki E, Yanagi S, Gotohda N, Kinoshita T, Kuwata T, Nagino M, Ochiai A, Fujii S: Intraoperative peritoneal lavage cytology offers prognostic significance for gastric cancer patients with curative resection.Cancer Sci 2017, 108:978-986.
7. Trilling B, Cotte E, Vaudoyer D, Isaac S, Piaton E, Maurice C, Glehen O, Passot G: Intraperitoneal-Free Cancer Cells Represent a Major Prognostic Factor in Colorectal Peritoneal Carcinomatosis.Dis Colon Rectum 2016, 59:615-622.
8. Rau B, Brandl A, Piso P, Pelz J, Busch P, Demtroder C, Schule S, Schlitt HJ, Roitman M, Tepel J, et al: Peritoneal metastasis in gastric cancer: results from the German database.Gastric Cancer 2020, 23:11-22.
9. Hamamoto Y: Complications in advanced or recurrent gastric cancer patients with peritoneal metastasis during and after palliative systemic chemotherapy.Mol Clin Oncol 2015, 3:539-542.
10. Huang Z, Liu D, Chen X, He D, Yu P, Liu B, Wu B, Hu J, Song B: Deep Convolutional Neural Network Based on Computed Tomography Images for the Preoperative Diagnosis of Occult Peritoneal Metastasis in Advanced Gastric Cancer.Front Oncol 2020, 10:601869.
11. Yang HK, Ji J, Han SU, Terashima M, Li G, Kim HH, Law S, Shabbir A, Song KY, Hyung WJ, et al: Extensive peritoneal lavage with saline after curative gastrectomy for gastric cancer (EXPEL): a multicentre randomised controlled trial.Lancet Gastroenterol Hepatol 2021, 6:120-127.
12. Zhang S, Zang D, Cheng Y, Li Z, Yang B, Guo T, Liu Y, Qu X, Che X: Identification of Key Gene and Pathways for the Prediction of Peritoneal Metastasis of Gastric Cancer by Co-expression Analysis.J Cancer 2020, 11:3041-3051.
13. Hinshaw DC, Shevde LA: The Tumor Microenvironment Innately Modulates Cancer Progression.Cancer Res 2019, 79:4557-4566.
14. Cristescu R, Lee J, Nebozhyn M, Kim KM, Ting JC, Wong SS, Liu J, Yue YG, Wang J, Yu K, et al: Molecular analysis of gastric cancer identifies subtypes associated with distinct clinical outcomes.Nat Med 2015, 21:449-456.
15. Barrett T, Troup DB, Wilhite SE, Ledoux P, Rudnev D, Evangelista C, Kim IF, Soboleva A, Tomashevsky M, Edgar R: NCBI GEO: mining tens of millions of expression profiles--database and tools update.Nucleic Acids Res 2007, 35:D760-765.
16. Davis S, Meltzer PS: GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor.Bioinformatics 2007, 23:1846-1847.
17. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.Proc Natl Acad Sci U S A 2005, 102:15545-15550.
18. Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS 2012, 16:284-287.
19. Hanzelmann S, Castelo R, Guinney J: GSVA: gene set variation analysis for microarray and RNA-seq data.BMC Bioinformatics 2013, 14:7.
20. Wilkerson MD, Hayes DN: ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking.Bioinformatics 2010, 26:1572-1573.
21. Team RC: R: A language and environment for statistical computing. 2013.
22. Oh SC, Sohn BH, Cheong JH, Kim SB, Lee JE, Park KC, Lee SH, Park JL, Park YY, Lee HS, et al: Clinical and genomic landscape of gastric cancer with a mesenchymal phenotype.Nat Commun 2018, 9:1777.
23. Tanaka Y, Chiwaki F, Kojima S, Kawazu M, Komatsu M, Ueno T, Inoue S, Sekine S, Matsusaki K, Matsushita H, et al: Multi-omic profiling of peritoneal metastases in gastric cancer identifies molecular subtypes and therapeutic vulnerabilities.Nature Cancer 2021, 2:962-977.
24. Wei C, Chen M, Deng W, Bie L, Ma Y, Zhang C, Liu K, Shen W, Wang S, Yang C, et al: Characterization of gastric cancer stem-like molecular features, immune and pharmacogenomic landscapes.Brief Bioinform 2021.
25. Mao D, Xu R, Chen H, Chen X, Li D, Song S, He Y, Wei Z, Zhang C: Cross-Talk of Focal Adhesion-Related Gene Defines Prognosis and the Immune Microenvironment in Gastric Cancer.Front Cell Dev Biol 2021, 9:716461.
26. Park HS, Kim YM, Kim S, Lee WS, Kong SJ, Yang H, Kang B, Cheon J, Shin SJ, Kim C, Chon HJ: High endothelial venule is a surrogate biomarker for T-cell inflamed tumor microenvironment and prognosis in gastric cancer.J Immunother Cancer 2021, 9.
27. Liu X, Yu H, Yan C, Mei Y, Lin C, Hong Y, Lin X, Zhang Q, Yu J: Plasmacytoid Dendritic Cells and ICOS(+) Regulatory T Cells Predict Poor Prognosis in Gastric Cancer: A Pilot Study.J Cancer 2019, 10:6711-6715.
28. Tucker MD, Brown LC, Chen YW, Kao C, Hirshman N, Kinsey EN, Ancell KK, Beckermann KE, Davis NB, McAlister R, et al: Association of baseline neutrophil-to-eosinophil ratio with response to nivolumab plus ipilimumab in patients with metastatic renal cell carcinoma.Biomark Res 2021, 9:80.
29. Zhang R, Li F, Li H, Yu J, Ren X: The clinical significance of memory T cells and its subsets in gastric cancer.Clin Transl Oncol 2014, 16:257-265.
30. Veen LM, Skrabanja TLP, Derks S, de Gruijl TD, Bijlsma MF, van Laarhoven HWM: The role of transforming growth factor beta in upper gastrointestinal cancers: A systematic review.Cancer Treat Rev 2021, 100:102285.
31. Koni M, Pinnaro V, Brizzi MF: The Wnt Signalling Pathway: A Tailored Target in Cancer.Int J Mol Sci 2020, 21.
32. Gong S, Duan Y, Wu C, Osterhoff G, Schopow N, Kallendrusch S: A Human Pan-Cancer System Analysis of Procollagen-Lysine, 2-Oxoglutarate 5-Dioxygenase 3 (PLOD3).Int J Mol Sci 2021, 22.
33. Low TY, van Heesch S, van den Toorn H, Giansanti P, Cristobal A, Toonen P, Schafer S, Hubner N, van Breukelen B, Mohammed S, et al: Quantitative and qualitative proteome characteristics extracted from in-depth integrated genomics and proteomics analysis.Cell Rep 2013, 5:1469-1478.

Tables

Table 1 Univariate COX regression based on DEGs