- 1Department of Surgical Oncology, The First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China
- 2Department of Medical Oncology, The First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China
- 3Department of Hematology & Oncology, The People's Hospital of Beilun District, Beilun Branch Hospital of the First Affiliated Hospital of Medical School of Zhejiang University, Ningbo, China
Background: Gastric cancer (GC) is the fifth most frequently diagnosed malignancy, and the third leading cause of tumor-related mortalities worldwide. Due to a high heterogeneity in GC, its treatment and prognosis are challenging, necessitating urgent identification of novel prognostic predictors for GC patients.
Methods: We downloaded RNA sequence data, from the Cancer Genome Atlas and microarray data from Gene Expression Omnibus database, then identified common differentially-expressed genes (DEGs) between GC and normal gastric tissues across four datasets. We then used a combination of protein-protein interaction (PPI) network and weighted gene co-expression network analysis (WGCNA) to identify key genes with prognostic value in GC. Thereafter, we used quantitative real time polymerase chain reaction (qRT-PCR) to validate expression of the identified key genes in the Zhejiang University (ZJU) cohort. Finally, we evaluated the relationships between gene expression and immune factors, including immune cells and biomarkers of immunotherapy.
Results: Among 426 common DEGs screened, 333 and 93 were upregulated and downregulated, respectively. PPI network and WGCNA successfully identified the top 30 hub genes, among which PTPRC, TYROBP, CCR1, CYBB, LCP2, and C1QB were common. Furthermore, TYROBP and C1QB were negatively associated with prognosis of GC patients, implying that they were key GC predictors. Interestingly, TYROBP and C1QB were positively correlated with predictive biomarkers for GC immunotherapy, including PD-L1 expression, CD8+ T cells infiltration, and EBV status.
Conclusions: TYROBP and C1QB were identified as two novel key genes with prognostic value in GC by network analysis.
Introduction
Gastric cancer (GC) is a major human health burden. According to the GLOBOCAN, in 2018 alone, there were over 1,000,000 new GC cases and an estimated 783,000 GC-related fatalities, making it the fifth prevalent cancer and the third leading cause of tumor mortality (1). Due to its high recurrence after surgery (2) and low sensitivity to chemotherapy (3), the overall 5-year survival rate of GC patients remains low. Therefore, it is urgent and crucial to identify novel prognostic biomarkers for GC patients.
With the rapid development and extensive application of high-throughput technology, vast amounts of gene expression profiles have been produced and utilized to identify differentially expressed genes (DEGs) by comparing tumor cells with the adjacent mucosa (4). However, the previous conventional studies have focused more on the individual DEGs while ignoring the complex network with a high degree of interconnection between the DEGs. Protein–protein interaction (PPI) networks and weighted gene co-expression network analysis (WGCNA) based on the microarray and RNA sequencing data have been shown to constitute powerful systematical biology strategies for mining the functional gene modules and identifying hub genes as candidate biomarkers, as well as therapeutic targets (5, 6). Over the past years, PPI and WGCNA have been extensively applied to screen out hub genes in multiple cancers. For instance, Chen et al. identified and validated that VCAN is associated with the progression and prognosis of pancreatic cancer by constructing a PPI network (7). Similarly, Yin et al. identified three novel blood-based diagnostic biomarkers for human hepatocellular carcinoma by WGCNA (8).
Herein, we constructed PPI and WGCNA networks based on the common DEGs from the TCGA-STAD (9) and 3 Gene Expression Omnibus (GEO) datasets [GSE65801 (10), GSE54129, and GSE118916 (11)]. Hub modules and hub genes were screened from the networks. An integrated bioinformatics analysis was performed to evaluate the function, pathway, and interrelation of the hub modules and the hub genes. We identified the key genes via survival analysis from the common hub genes derived from the PPI and WGCNA network, then validated them in the Oncomine database, ZJU cohort, and GSE15459 dataset (12). An immune analysis was performed to investigate the association between the key genes and the immune factors using the TCGA-STAD and GSE51575 dataset (13).
Materials and Methods
Study Design
The design of this study is shown as a workflow (Figure 1). We screened the differentially expressed genes (DEGs) between the GC and normal or adjacent mucosa tissue from the four cohort profile datasets, i.e., TCGA-STAD (9), GSE65801 (10), GSE54129, and GSE118916 (11). The construction of the protein–protein interaction (PPI) network and weighted gene co-expression network was based on the DEGs, and we identified the common hub genes from the networks. The expression of the common hub genes was validated in the Oncomine database and ZJU cohort. We performed survival analyses of the common hub genes using the TCGA-STAD dataset. Immune analyses were performed to evaluate the correlation between the key genes and the tumor microenvironment using the TCGA-STAD and GSE51575 datasets (13).
Figure 1. Workflow of our study for identifying key genes with prognostic value in gastric cancer, including data preparing, processing and analysis.
Data Collection
We downloaded the RNA sequencing data and clinical datasets of GC patients from the TCGA repository of the National Cancer Institute (https://cancergenome.nih.gov/). The TCGA-STAD datasets constituted 375 tumor and 32 normal samples. Microarray data of GSE65801, GSE54129, GSE118916, GSE15459, and GSE51575 datasets were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). The GSE65801 microarray data was downloaded from the GPL14550 Platform (Agilent-028004 SurePrint G3 Human GE 8x60K Microarray, Probe Name version, Agilent Technologies) and included 32 gastric cancer tissues and 32 paired noncancerous tissues (Submission date: Feb 10, 2015) (10). The microarray data of GSE54129 and GSE15459 was downloaded from the GPL570 Platform ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array), with the former constituting 111 gastric cancer tissues and 21 normal gastric tissues (Submission date: Jan 16, 2014) while the latter included 200 primary gastric cancer tissues (Submission date: Mar 30, 2009) (12). The GSE118916 microarray data was downloaded from the GPL15207 Platform ([PrimeView] Affymetrix Human Gene Expression Array) and included 15 gastric cancer tissues and 15 paired adjacent mucosa tissues (Submission date: Aug 22, 2018) (11). The GSE51575 microarray data was downloaded from the GPL13607 Platform (Agilent-028004 SurePrint G3 Human GE 8x60K Microarray, Feature Number version) and included 26 adjacent mucosa tissues, 14 EBV-positive gastric cancer tissues, and 12 EBV-negative gastric cancer tissues (Submission date: Oct 23, 2013) (13). The GSE51575 dataset was derived from a primary study (13) and contained some essential information for our research, including gene expression of immune checkpoints, and EBV infection status. The acquisition and application methods of all the data were according to the guidelines and policies of the GEO and TCGA databases.
Data Preprocessing and Common DEGs Identification
The retrieved gene expression data from the GEO database was preprocessed, including background correction and normalization in the R version 3.6.1 software. We utilized the Bioconductor Annotation Data software package to transform the microarray data probes to gene symbols. When several probes were matched to the same gene symbol, the median value was set as the final expression value of the gene. The “limma” and “edgeR” R packages were utilized to identify the DEGs between the GC tissues and normal or adjacent mucosa tissues in the GEO and TCGA datasets, respectively (14, 15). Genes with adjusted P < 0.05 and |Fold change (FC)| > 1.5 were selected as the DEGs. Common DEGs were defined as the overlap of the DEGs from the TCGA-STAD, GSE65801, GSE54129, and GSE118916. The Venn diagram was generated online (http://bioinformatics.~psb.ugent.be/webtools/Venn/).
Functional Annotation, Pathway Enrichment, and Interrelation Analysis
We analyzed the functional annotation and pathway enrichment using the Database for Annotation, Visualization and Integrated Discovery (DAVID) web portal (https://david.ncifcrf.gov/) (16, 17). After uploading the list of common DEGs, we obtained the Gene Ontology (GO) enrichment results of the biological process (BP), cellular component (CC), molecular function (MF), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Interrelation analysis between pathways and hub genes was performed using the ClueGo (Version 2.5.4), a plug-in of Cytoscape software (18). P < 0.05 was set as the cut-off criteria.
PPI Network Construction and MCODE Analysis
First, we utilized the STRING database (http://string-db.org) to construct the PPI network of DEGs and interactions, with a combined score > 0.4 considered statistically significant (6). After that, we used the Cytoscape software (Version 3.7.2) to visualize the PPI network (19). Subsequently, the Molecular Complex Detection (MCODE, version 1.5.1) plug-in tool of Cytoscape was used to screen and visualize the hub modules in the PPI network with the MCODE score = 5, degree = 2, Node score cut-off = 0.2, K-score = 2, and Max. Depth = 100 (20). The functional annotation for the genes in the modules was performed using the DAVID.
WGCNA Network Construction
We utilized the WGCNA to analyze the co-expressed gene module and identify the hub module correlated to the clinical traits (5). In this study, we selected the common DEGs for the WGCNA network construction using the “WGCNA” R package. Sample clustering of the common DEGs was applied to filter the outlier sample with a height cut-off value of 20,000. A power of β = 4 and minimum module size = 30 were set as per the standard scale-free networks. The adjacencies between all the filtered genes were conducted and converted into a topological overlap matrix (TOM) and the corresponding dissimilarity (1-TOM). The hierarchical clustering function was used to classify the genes with a high absolute correlation into modules based on the TOM-based dissimilarity for the gene dendrogram. The dissimilarity of the module eigengenes was calculated to merge similar modules with a height cut-off value of 0.25. Module eigengene (ME), defined as the first principal component of a given module, was regarded as the representative of the module. The correlation between the ME and the clinical traits, including age, gender, grade, and the stage, was one of the factors for identifying the hub module. Gene significance (GS) was defined as the log10 transformation of P-value in the linear regression between the gene expression and the clinical traits. The module membership (MM) was identified as the correlation between the gene expression and the ME. The hub module was identified by the highest correlation between the ME and the clinical traits, as well as the most significant correlation between the MM and the GS. Subsequently, the hub module was visualized using the Cytoscape software. The functional annotation and pathway enrichment for the genes in the hub module were conducted using DAVID.
Common Hub Genes Identification and Validation
The highly interconnected hub genes with the other nodes in a module were regarded as functionally significant genes. Herein, the hub modules were identified using the PPI network and the WGCNA network. The hub genes in the hub modules were screened using the cytoHubba (Version 0.1) tool, a plug-in of the Cytoscape software (21). The hub genes that ranked the top 30 in the hub modules were selected as the candidates using the Degree method, and the interrelation analysis was performed, as described previously. The common hub genes defined as the overlap of the hub genes from the PPI network and the WGCNA network were identified for further analysis and validation. The expression of the common hub genes was validated using the Oncomine database.
Tissue Samples and Total RNA Isolation
We obtained 10 pairs of the GC tissues and adjacent mucosa tissues from GC patients who underwent surgery at the First Affiliated Hospital of Zhejiang University (ZJU cohort), excluding those who had been exposed to pre-operative chemotherapy or radiotherapy. The Institutional Review Board of the First Affiliated Hospital of Zhejiang University approved the protocol of this study. All the GC patients signed informed consent. The total RNA from each of the 10 GC tissues and 10 paired adjacent mucosa tissues was isolated using a RNeasy Mini Kit (Cat.no.74106, Qiagen, Germany) and quantified using a NanoDrop One (Cat. ND-ONE-W, ThermoFisher Scientific, USA).
Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Briefly, we performed reverse transcription to synthesize the first-strand cDNA using 1 μg total RNA isolated from the GC tissues and paired normal mucosa tissue samples using the PrimeScript™ RT Master Mix (Perfect Real Time) (Cat. #RR036A, TaKaRa, Japan). After that, qRT-PCR was performed using the TBGreen®Premix Ex Taq™II (Tli RNase H Plus) (Cat. #RR820A, TaKaRa, Japan). We utilized the Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene as an internal control.
Survival Analysis and Key Genes Identification
We assessed the prognostic value of the common hub genes by evaluating the association of the gene expression and overall survival of patients in the TCGA-STAD dataset. The key genes with clinical significance were identified through expression validation and survival analysis of the common hub genes. Finally, the transcriptional expression of the key genes was validated in the ZJU cohort.
Immune Analysis of Key Genes
As an interactive web platform, Tumor Immune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/) is utilized to estimate tumor immune infiltration across diverse cancer types (22). In this study, we analyzed the correlation of the expression of key genes in GC with immune infiltration (CD8+ T Cells, CD4+ T Cells, and Macrophages) using the “Gene” module and explored the correlation between immune infiltration and survival using the “Survival” module in TIMER. CIBERSORT is a deconvolution algorithm used to estimate the immune cell type proportions with a signature matrix of 547 genes by support vector regression (23). The output includes a P-value for the deconvolution of each sample using the Monte Carlo sampling after running with 1,000 permutations. The CIBERSORT P-value reflects the statistical significance of the results, and a threshold < 0.05 is recommended. We uploaded the gene expression profile constituting 375 tumor samples in the TCGA-STAD to the CIBERSORT web portal (https://cibersort.stanford.edu/). Consequently, 240 samples with CIBERSORT P < 0.05 were included in calculating the Spearman's correlation between the key genes and 22 types of the infiltrating immune cells. Besides, we evaluated the correlation between the key genes and the immune factors, including cytolytic activity molecules (GZMA and PRF1) and the immune checkpoints (CD274, PDCD1, PDCD1LG2, VTCN, and LAG3) of 375 samples in the TCGA-STAD dataset. The correlation between the key genes and the immune checkpoints was validated using the GSE51575 data. Information about the TCGA-STAD subtypes, including CIN, EBV, MSI, and GS, was mined from the cBioportal database (24), an open-access resource providing data from the TCGA project (https://www.cbioportal.org/). Epstein Barr virus (EBV) associated with gastric cancer was classified as one of the four molecular subtypes in 2014 (9). EBV positive status additionally was one of the biomarkers for immunotherapy. We then evaluated the association between the key genes and the EBV status in the TCGA-STAD and GSE51575 datasets. The correlation between the key genes and immune factors in the TCGA-STAD dataset was visualized using the MeV (MultiExperiment Viewer, version 4.9.0.) software (25).
Correlation Analysis of Key Genes
The correlation analysis of the key genes was performed using the Gene Expression Profiling Interactive Analysis (GEPIA) (26), a web server for analyzing gene expression from the TCGA and GTEx samples (http://gepia2.cancer-pku.cn/).
Statistical Analysis
Statistical analyses were performed using the SPSS 21.0 software. We used the Kaplan-Meier survival to analyze the association between gene expression and the overall survival. The log-rank test was used to determine significant differences in the survival curves stratified by the gene expression level. We calculated the median overall survival time, and the 95% confidence interval where relevant. The correlations of the gene expression with immune cells and immune factors were evaluated using the Spearman's correlation and statistical significance. The continuous variables in the two groups and multi-subgroups were compared using the Student's t-test and ANOVA, respectively. A P < 0.05 was considered statistically significant.
Results
Identification, Functional Annotation, and Pathway Analysis of Common DEGs
The common DEGs between the GC tumor and normal or adjacent mucosa tissues were screened from the GEO datasets (GSE65801, GSE54129, and GSE118916), and the TCGA-STAD dataset. Consequently, 426 common DEGs were identified, of which 333 were upregulated, and 93 were downregulated (Figures 2A,B; Table S1). After that, we performed Go analyses by uploading the identified 426 common DEGs to the DAVID web portal. The DEGs were enriched as per the four subontologies: BP, CC, MF, and the KEGG pathway. For BP (Figure 2C), the DEGs were primarily enriched in signal transduction, cell adhesion, and immune response. The other biological processes were associated with the tumor microenvironment constituted the extracellular matrix organization, collagen catabolic process, collagen fibril organization, inflammatory response, leukocyte migration, angiogenesis, chemotaxis, and cell-cell signaling. For MF (Figure 2D), enrichment of the DEGs was primarily in protein binding and calcium ion binding. The other molecular functions were mainly enriched in receptor activity included integrin binding, heparin binding, collagen binding, platelet-derived growth factor binding, cytokine activity, identical protein binding, extracellular matrix binding, protease binding, glycosaminoglycan binding, indanol dehydrogenase activity, and the N-formyl peptide receptor activity. For CC (Figure 2E), the DEGs were mainly involved in the plasma membrane, extracellular space, and exosome. The other cellular components included extracellular region, extracellular matrix, proteinaceous extracellular matrix, collagen trimer, integral component of the plasma membrane, cell surface, endoplasmic reticulum lumen, plasma membrane, external side of the plasma membrane, basement membrane, membrane, membrane raft, podosome, and the cell projection. For the KEGG pathway (Figure 2F), the DEGs were mainly enriched in the PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and focal adhesion. The other pathways were mainly enriched in immune signaling, including the TNF signaling pathway, staphylococcus aureus infection, rheumatoid arthritis, leukocyte transendothelial migration, and chemical carcinogenesis.
Figure 2. Identification, functional annotation and pathway analysis of the common DEGs in four cohorts (TCGA-STAD, GSE65801, GSE54129, and GSE118916). (A) Venn diagram of the up-regulated genes in the four cohorts. (B) Venn diagram of the down-regulated genes in the four cohorts. (C) Biological processes of the common DEGs. (D) Cellular components of the common DEGs. (E) Molecular functions of the common DEGs. (F) KEGG pathway analysis of the common DEGs.
Identification, Enrichment, and Interrelation Analysis of Hub Modules and Hub Genes in the PPI Network
The PPI network was constructed using the STRING online database, and the top 3 significant modules were screened using the MCODE (Figure 3A). Module 1 contained 36 genes, with 35 upregulated and 1 down regulated in the tumor tissues vs. the normal tissues. Module 2 consisted of 38 genes, which were upregulated in the tumor tissues. Module 3 constituted 28 genes, with 27 upregulated and 1 downregulated in the tumor tissues. We performed function annotation of the modules using DAVID, as described previously. The GO analysis results (Figure 3B; Table S2) disclosed that module 1 was primarily enriched in the extracellular matrix organization, extracellular matrix structural constituent, and the extracellular region. Module 2, on the other hand, was primarily enriched in signal transduction, protein binding, and extracellular space, whereas module 3 was mainly involved in the inflammatory response, protein binding, and the plasma membrane. In addition, the top 30 hub genes with a high degree of connectivity in the PPI network were identified using cytoHubba (Supplementary Figure 1; Tables S3, S4). We performed an interrelation analysis between the pathways in the BPs of the hub genes using ClueGo to evaluate the pathway enrichment of the hub genes and the crosstalk between pathways. Consequently, the hub genes were primarily enriched in positive regulation of the response to macrophage colony-stimulating factor, positive regulation of the tumor necrosis factor biosynthetic process, negative regulation of the myeloid cell apoptotic process, and fibrillar collagen trimer (Figure 3C; Supplementary Figure 2). Based on the results above, we observed an enrichment of the hub modules and hub genes in the inflammatory response and extracellular matrix organization.
Figure 3. Identification, enrichment and interrelation analysis of the hub modules and the hub genes in PPI network. (A) Identification of the top 3 hub modules in the PPI network. (B) GO term enrichment analysis of the hub modules. (C) Interrelation analysis in the biological process pathways of the hub genes.
Identification of the Hub Module and Hub Genes in WGCNA Network
We constructed the WGCNA network using the “WGCNA” R package. The expression patterns of the genes in the same module were similar and relevant by the average linkage clustering. We included 315 samples with clinical traits to filter the outlier samples via sample clustering of the common DEGs, and 17 samples were excluded with the height of 20,000 (Supplementary Figure 3A). A soft threshold (β) = 4 was set to ensure a scale-free network (R2 = 0.94; Supplementary Figures 3B–E). Similar modules with a height cut-off value of 0.25 were merged (Figure 4A), and 3 modules marked in blue, turquoise, and gray were identified (Figure 4B). The blue module contained 112 genes and the turquoise module 188 genes. Besides, 126 genes not included in any module were put into the gray module. The gray module was identified as not co-expressed and would be excluded in subsequent analyses. The interaction of the modules was visualized as the network heatmap (Figure 4C), which indicated that genes in the same module had a highly co-expressed relationship with each other. Then, the correlation between the GS and MM in the turquoise and blue modules was calculated, respectively. The correlation was significant in the blue module (R = 0.81, P < 0.001; Figure 4D) and not in the turquoise module (R = 0.14, P = 0.054; Figure 4E). Furthermore, the relationship between the modules and the clinical traits was evaluated to identify the hub module. The result showed that the blue module was significantly associated with the GC grade (R = 0.31, P < 0.001; Figure 4F). The top 30 hub genes in the blue module were screened using cytoHubba via the Degree method (Supplementary Figure 4; Table S5, Table S6). Consequently, the blue module was identified as the hub module in the WGCNA network and the top 30 hub genes in the blue module.
Figure 4. Construction and identification of the WGCNA co-expression modules associated with the clinical traits, based on the common DEGs expression data of TCGA-STAD. (A) Cluster dendrogram of the module eigengenes. The dissimilarity of module eigengenes is calculated to merge some similar modules with a height cut-off value of 0.25. (B) Cluster dendrogram of the DEGs. Highly similar modules are identified by clustering and then merged dynamically. (C) Network heatmap plot of the common DEGs. The branch in the hierarchical clustering dendrograms corresponds to each module. Color bars beneath and toward the right of the dendrograms show the color-coded module membership. The more saturated yellow and red indicate the higher co-expression interconnectedness in the heatmap. (D) Scatter plot of the GS for the grade vs. the MM in the blue module. (E) Scatter plot of the GS for the grade vs. the MM in the turquoise module. (F) Heatmap of the correlation between the module eigengenes and the clinical traits of gastric cancer. The blue module is the most positively correlated with the grade and identified as the hub module. GS, gene significance; MM, module membership.
Function Annotation, Pathway Enrichment, and Interrelation Analysis of the Blue Module and the Hub Genes in the WGCNA Network
Function annotation and pathway enrichment of the blue module in the WGCNA network were performed in DAVID, as previously described. It was mainly enriched in the immune response and cell adhesion for BP, plasma membrane for CC, and protein binding for MF (Figure 5A). The KEGG pathway enrichment analysis identified the cytokine-cytokine receptor interaction as the most significantly enriched pathway, and the other pathways included the chemokine signaling pathway, tuberculosis, phagosome, and osteoclast differentiation (Figure 5B). Furthermore, an interrelation analysis between the pathways in the BPs of the hub genes was performed, as previously described. The hub genes were primarily involved in positive macrophage activation, regulation of tumor necrosis factor production, and regulation of tolerance induction (Figure 5C; Supplementary Figures 5A,B).
Figure 5. Function annotation and pathway enrichment of the hub module and hub genes in the WGCNA network. (A) GO analysis of the hub module in the WGCNA network. (B) KEGG pathway analysis of the hub module in the WGCNA network. (C) Interrelation analysis in the biological process pathways of the hub genes.
Identification, Validation, and Survival Analysis of the Key Genes
The top 30 hub genes in the PPI and WGCNA networks were overlapped to identify the common hub genes, including PTPRC, TYROBP, CCR1, C1QB, CYBB, and LCP2 (Figure 6A). We then conducted a literature review to investigate the association between the hub genes in the networks and the tumors. Consequently, among the common hub genes, 50.0% (3/6) genes have been shown to promote tumor progression in gastric cancer. Among the other hub genes in the PPI and WGCNA networks, 70.8% (17/24) and 29.2% (7/24) have been reported in gastric cancer-associated studies, respectively (Table S7). We next focused on the common hub genes that might play a vital role in gastric cancer, considering the strong connection between the common hub genes and the other hub genes in the PPI and WGCNA networks. Based on the Oncomine database, using Data type = mRNA, P-value < 0.05, |FC| > 1.5 and gene rank = “all” as the threshold, the expression levels of the common hub genes were significantly higher in the GC tumor tissues compared with the normal or adjacent mucosa tissues (Figure 6B). Then, survival analyses of the common hub genes were performed using the TCGA-STAD dataset. The results showed that TRYOBP and C1QB were negatively associated with the overall survival of the GC patients (PTYROBP = 0.029, PC1QB = 0.030; Figure 6C), which was validated using the GSE15459 dataset (PTYROBP = 0.001, PC1QB = 0.001; Supplementary Figures 6A,B). However, CCR1, CYBB, LCP2, and PTPRC were not associated with the overall survival (PCCR1 = 0.412, PLCP2 = 0.148, PPTPRC = 0.132, PCYBB = 0.189; Supplementary Figures 7A–D). Based on these results, we further identified TYROBP and C1QB as the two key genes with prognostic value in gastric cancer. In addition, the high expression of TYROBP and C1QB were also validated in the GC tumor tissues compared with the adjacent mucosa tissues using the ZJU cohort (PTYROBP = 0.045, PC1QB = 0.031; Figure 6D). Besides, we explored the expression location of TYROBP and C1QB using the Human Protein Atlas database. As shown in the Supplementary Figure 8, TYROBP and C1QB were mainly located in tumor cells, which needed further validation in the experiments. Furthermore, we investigated the prognostic value of the other hub genes (the PPI or WGCNA network contained 24 genes) via univariate Cox analysis. We identified 50% (12/24) of the hub genes in the PPI network and 16.7% (4/24) of the hub genes in the WGCNA network that were negatively associated with the overall survival in GC patients (Table S8).
Figure 6. Identification, validation and survival analysis of key genes. (A) Identification of the common hub genes between the PPI network and the WGCNA network. (B) Validation of the common hub genes in ONCOMINE database. Red color represents relatively higher expression of the common hub genes in tumors than normal tissues. Numbers represented the number of studies. (C) Survival analysis of the key genes in TCGA-STAD dataset. (D) Validation of the key genes at the transcriptional level in ZJU cohort. *P < 0.05.
The Correlation Between TYROBP, C1QB, and Immune Factors, Including Immune Cells and Biomarkers for Immunotherapy
As mentioned, we found that TYROBP and C1QB could be involved in the immune response via the function annotation, pathway enrichment, and interrelation analysis. Therefore, we evaluated the correlation between TYROBP, C1QB, and immune factors, respectively. The TIMER database analysis results revealed a significantly positive correlation of both genes with the CD8+ T cells, CD4+ T cells, and macrophages (Figure 7A). Survival analyses of the immune cells were also performed using the TIMER database, and the results indicated that macrophages were negatively associated with the survival time of GC patients (P = 0.004; Figure 7B). Furthermore, using the CIBERSORT algorithm, we estimated the proportion of the infiltrating immune cells in the tumor microenvironment. Consequently, we found that the infiltrating immune cells mainly consisted of the CD8+ T cells, CD4+ T cells, and macrophages. TYROBP and C1QB were negatively correlated with resting memory CD4+ T cells (R1 = −0.27, R2 = −0.25; Figure 7C; Table S9) and positively correlated with the CD8+ T cells (R1 = 0.21, R2 = 0.26; Figure 7C; Table S9), activated memory CD4+ T cells (R1 = 0.26, R2 = 0.37; Figure 7C; Table S9), and macrophage M2 (R1 = 0.46, R2 = 0.47; Figure 7C; Table S9). We performed correlation analyses for TYROBP and C1QB with differential markers of macrophages to further investigate the association of TYROBP and C1QB with the macrophages. As shown in the Supplementary Figure 9, the correlation of CD11b and CD206 (M2) was stronger compared with the CD68 (M0) and CCR7 (M1). Furthermore, we performed univariate and multivariate Cox regression for TYROBP, C1QB, and macrophages, respectively. As shown in Table S10, TYROBP (HR = 1.455, P = 0.029), C1QB (HR = 1.474, P = 0.030), and macrophage M2 (HR = 1.494, P = 0.024) were identified as the significant risk factors, while macrophage M1 (HR = 1.182, P = 0.374) was not associated with the overall survival via the univariate Cox regression analysis. Notably, in the multivariate Cox regression model, both TYROBP (HR = 1.399, P = 0.073) and C1QB (HR = 1.428, P = 0.067) were not significantly correlated with the overall survival when adjusted by macrophage M2. However, TYROBP (HR = 1.518, P = 0.023) and C1QB (HR = 1.550, P = 0.024) were still significantly associated with poor outcomes when adjusted by macrophage M1. These findings suggested that macrophage M2 is involved in TYROBP/C1QB-mediated progression and poor survival outcomes in GC. Then, the cytolytic activity of the immune cells was estimated using the average expression of GZMA and PRF1 (27). The results showed that TYROBP and C1QB were negatively correlated with the cytolytic activity (R1 = −0.30, R2 = −0.28; Figure 7D; Table S9), which revealed the immunosuppressive microenvironment in tumors. Furthermore, the correlation between TYROBP, C1QB, and immune checkpoints was assessed, respectively. Consequently, TYROBP and C1QB were positively correlated with CD274, PDCD1, PDCD1LG2, and LAG3 but negatively correlated with VTCN1 in the TCGA-STAD dataset (Figure 7E; Table S9) and GSE51575 dataset (Figure 7F). In addition, we evaluated the expression of TYROBP and C1QB in the TCGA subtypes. Compared with the other subtypes, the expression of TYROBP and C1QB in the EBV positive subtype was significantly higher (PTYROBP < 0.0001, PC1QB < 0.0001; Figure 7G). These results were validated using the GSE51575 dataset (PTYROBP < 0.0001, PC1QB = 0.0002; Figure 7H). In summary, we established that the expression of TYROBP and C1QB was positively correlated with the PD-L1 expression, CD8+ T cell infiltration, and the EBV status; three predictive biomarkers for immunotherapy.
Figure 7. Correlation analysis between TYROBP, C1QB, and immune factors, including immune cells and biomarkers for immunotherapy. (A) Correlation between TYROBP, C1QB, and the immune cells in the TIMER database. (B) Survival analysis of the immune cells in the TIMER database. (C) Correlation between TYROBP, C1QB, and the immune cells in the TCGA-STAD dataset. (D) Correlation between TYROBP, C1QB, and the immune cytolytic activity in the TCGA-STAD dataset. (E) Correlation between TYROBP, C1QB, and the immune checkpoints in the TCGA-STAD dataset. (F) Correlation between TYROBP, C1QB and the immune checkpoints in the GSE51575 dataset. (G) Relative expression of TYROBP and C1QB in the four molecular subtypes in the TCGA-STAD dataset. (H) Relative expression of TYROBP and C1QB in normal tissues, EBV negative (–) and positive (+) GC tissues in the GSE51575 dataset. R1, the correlation coefficient between TYROBP and the immune factors; R2, the correlation coefficient between C1QB and the immune factors; GC, gastric cancer. The data marked red color is statistically significant. ***P < 0.001, ****P < 0.0001.
Discussion
GC is the third leading cause of global cancer-related deaths. However, to date, effective treatments have not yet been developed, owing to a limited understanding of the molecular mechanisms underlying GC development. Over the past years, the applications of PD-1/PD-L1 checkpoint blockades in cancer have revolutionized oncology (28, 29). Particularly, these approaches have guided immunotherapy strategies against multiple cancers, such as melanoma (30), lung cancer (31), glioblastoma (32) and liver cancer (33). However, the efficacy and responsiveness of immunotherapeutic agents significantly varies across GC patients, largely due to high tumor heterogeneity and molecular complexity (34). Thus, it is crucial to unravel the underlying molecular mechanisms of GC tumorigenesis and progression and identify potential prognostic and therapeutic targets. In the present study, analysis of four gene expression profiles revealed common DEGs, mainly enriched in signal transduction, cell adhesion and immune response. Previous studies have shown that extracellular matrix remodeling and abnormal immune microenvironment play important roles in tumorigenesis and tumor progression. In the past decade, studies have demonstrated that the interaction between tumor microenvironment and tumor cells is essential for tumor biological behavior (35–38). In the present study, PPI and WGCNA networks revealed 6 common hub genes, including PTPRC, TYROBP, CCR1, C1QB, CYBB, and LCP2. According to Wang et al. (39), CYBB is associated with invasion and prognosis of human gastric cancer, whereas PTPRC, also known as CD45, has been previously used to assess the extent of immune cell infiltration in intestinal-type Japanese gastric cancer (40). On the other hand, Chen et al. (41) previously reported that CCR1 was associated with CD4+CD25+ Tregs of regional lymph nodes in forestomach carcinoma. Interestingly, TYROBP and C1QB were both correlated with immune infiltration levels, suggesting a potential key role in prognosis of GC patients. These factors have previously been positively associated with three predictive biomarkers for immunotherapy in GC, including PD-L1 expression (42, 43), CD8+ T cells infiltration (44) and EBV status (45, 46).
Previous studies have shown that TYROBP, also known as DAP12, is overexpressed and related to tumor progression in multiple cancers. Functionally, its encoded protein, a transmembrane signaling polypeptide on the surface of a variety of immune cells, mediates signaling transductions (47, 48). For example, Shabo et al. (49) reported an association between high TYROBP expression with skeletal and liver metastases as well as poor survival of breast cancer patients. Similarly, Cheray et al. (50) implicated TYROBP in glioblastoma tumorigenesis and aggressiveness. In the present study, TYROBP overexpression was associated with poor survival of GC patients. In addition, results from interrelation analysis showed that TYROBP was associated with positive macrophage activation, regulation of tumor necrosis factor production and regulation of tolerance induction. This is consistent with a previous study that found a positive association between TYROBP with macrophage M2, as well as the immunosuppressive and pro-tumorigenic subtype of macrophage in the tumor microenvironment (51). Similarly, Takamiya et al. (52) found that TYROBP was involved in the interaction between lung cancer cell and macrophage M2 to enhance TGF-β secretion in vitro. Our results further revealed a positive correlation between TYROBP and CD8+ T cells, but a negative association with cytolytic activity. In addition, we found a positive association between TYROBP with most checkpoints, including CD274, PDCD1, PDCD1LG2, and LAG3. Taken together, these results indicated that TYROBP might be playing an immunosuppressive role on CD8+ T cells and macrophages to promote tumor immune escape in gastric cancer. Coincidentally, Yoshida et al. (53) reported that TYROBP deficiency in liver allografts resulted in activation of graft-infiltrating CD8+ T cells and production of pro-inflammatory cytokine, whereas Kovats et al. (54) found that loss of TYROBP and FcRγ promoted IL-12 production and CD8+ T cell response by CCR2+ Mo-DCs. Thus, TYROBP might be a negative factor in anti-tumor immune response. Furthermore, we found a significantly higher TYROBP expression in EBV positive patients relative to those with other subtypes. To date, EBV status is one of the validated predictive biomarkers for immunotherapy (45, 46). These results suggest that TYROBP might be associated with the multiple biomarkers for immunotherapy in gastric cancer, although further validation using large clinical cohorts is required.
In the present study, we found significantly higher expression of C1QB, that encodes the C1qB chain, in tumor than adjacent normal GC tissues. In addition, C1QB was negatively associated with prognosis of GC patients. Previous studies have shown that C1q, the first recognition subcomponent of the complement classical pathway, comprises three chains (C1qA, C1qB, and C1qC) and exerts complex effects on tumorigenesis of multiple tumors, including prostate (55), and ovarian cancer (56) as well as gliomas (57). Yamada et al. (58) reported that high C1QB expression was significantly related to poor prognosis in renal cell carcinoma. On the other hand, Linnartz-Gerlach et al. (59) found that C1qB was downregulated in the brain of triggering receptor expressed on myeloid cells-2 (TREM2) knock-out mice. Interestingly, TREM2 has been reported to transmit intracellular signals through the associated transmembrane adapter TYROBP (60). In the present study, we found a strong correlation between TYROBP and C1QB expression in GC patients (R = 0.92, P < 0.001; Supplementary Figure 10). Studies have shown that a dysregulation of this signaling pathway leads to a wide range of pathophysiological changes and diseases, such as aging (59), bone cysts (61) and Alzheimer's disease (62). In our study, we also found an association between C1QB with PD-L1 expression, CD8+ T cells infiltration and EBV status, which was very similar to the TYROBP pattern. However, in vitro and in vivo studies are needed to validate the observed relationship between TYROBP and C1QB in GC patients.
This study had some limitations. Firstly, our results will be more convincing and interesting through additional validation of TYROBP and C1QB in vivo in vitro experiment. For example, immunofluorescent detection is more precise than immunohistochemistry for analyzing co-localization of TYROBP and C1QB. Secondly, although our integrated network analysis indicated the prognosis value of TYROBP and C1QB in gastric cancer, further validation is needed using more clinical cohorts. Importantly, our bioinformatics findings indicated that macrophage M2 might be involved in TYROBP/C1QB-mediated progression and poor survival outcomes in GC. Further experimental studies are needed to unravel the role of macrophage M2 in GC.
In conclusion, we used integrated network analysis, PPI and WGCNA, to reveal overexpression of TYROBP and C1QB, and affirm their prognostic value in GC patients. To our knowledge, this is the first report associating TYROBP and C1QB with GC progression and prognosis. Our findings lay a foundation for future research aiming to elucidate the role of these genes in GC tumorigenesis and progression.
Data Availability Statement
Publicly available datasets were analyzed in this study, these can be found in The Cancer Genome Atlas (https://portal.gdc.cancer.gov/); the NCBI Gene Expression Omnibus (GSE65801, GSE54129, GSE15459, GSE51575, and GSE118916).
Ethics Statement
The studies involving human participants were reviewed and approved by Research Ethics Committee of the First Affiliated Hospital, College of Medicine, Zhejiang University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
LT conceived and designed the study. JJ and YD collected data, performed data analysis, and wrote manuscript. MW, YC, HaifW, and HaiyW were involved in data interpretation and critically reviewed the manuscript. XL performed qRT-PCR. All authors read and approved the final manuscript.
Funding
This study was supported by Grants from the Key Program of Natural Science Foundation of Zhejiang (No. LZ16H160001), Grants from the Natural Science Foundation of Zhejiang (No. LQ20H160043), the Key Project of Scientific and Technological Innovation of Zhejiang Province (No. 2018C03022), and the Key Project jointly built by Zhejiang province and Health Planning Committee (No. 2017209495).
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.
Acknowledgments
All supports from participants in this research were undeniable.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.01765/full#supplementary-material
References
1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2018) 68:394–424. doi: 10.3322/caac.21492
2. Isobe Y, Nashimoto A, Akazawa K, Oda I, Hayashi K, Miyashiro I, et al. Gastric cancer treatment in Japan: 2008 annual report of the JGCA nationwide registry. Gastric Cancer. (2011) 14:301–16. doi: 10.1007/s10120-011-0085-6
3. Cunningham D, Stenning SP, Smyth EC, Okines AF, Allum WH, Rowley S, et al. Peri-operative chemotherapy with or without bevacizumab in operable oesophagogastric adenocarcinoma (UK Medical Research Council ST03): primary analysis results of a multicentre, open-label, randomised phase 2–3 trial. Lancet Oncol. (2017) 18:357–70. doi: 10.1016/S1470-2045(17)30043-8
4. Chatsirisupachai K, Palmer D, Ferreira S, de Magalhaes JP. A human tissue-specific transcriptomic analysis reveals a complex relationship between aging, cancer, and cellular senescence. Aging Cell. (2019) 18:e13041. doi: 10.1111/acel.13041
5. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559
6. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res. (2015) 43:D447–52. doi: 10.1093/nar/gku1003
7. Chen Q, Yu D, Zhao Y, Qiu J, Xie Y, Tao M. Screening and identification of hub genes in pancreatic cancer by integrated bioinformatics analysis. J Cell Biochem. (2019) 120:19496–508. doi: 10.1002/jcb.29253
8. Yin L, He N, Chen C, Zhang N, Lin Y, Xia Q. Identification of novel blood-based HCC-specific diagnostic biomarkers for human hepatocellular carcinoma. Artif Cells Nanomed Biotechnol. (2019) 47:1908–16. doi: 10.1080/21691401.2019.1613421
9. Cancer Genome Atlas Research N. Comprehensive molecular characterization of gastric adenocarcinoma. Nature. (2014) 513:202–9. doi: 10.1038/nature13480
10. Li H, Yu B, Li J, Su L, Yan M, Zhang J, et al. Characterization of differentially expressed genes involved in pathways associated with gastric cancer. PLoS ONE. (2015) 10:e0125013. doi: 10.1371/journal.pone.0125013
11. Li L, Zhu Z, Zhao Y, Zhang Q, Wu X, Miao B, et al. FN1, SPARC, and SERPINE1 are highly expressed and significantly related to a poor prognosis of gastric adenocarcinoma revealed by microarray and bioinformatics. Sci Rep. (2019) 9:7827. doi: 10.1038/s41598-019-43924-x
12. Ooi CH, Ivanova T, Wu J, Lee M, Tan IB, Tao J, et al. Oncogenic pathway combinations predict clinical prognosis in gastric cancer. PLoS Genet. (2009) 5:e1000676. doi: 10.1371/journal.pgen.1000676
13. Kim SY, Park C, Kim HJ, Park J, Hwang J, Kim JI, et al. Deregulation of immune response genes in patients with epstein-barr virusassociated gastric cancer and outcomes. Gastroenterology. (2015) 148:137–47 e9. doi: 10.1053/j.gastro.2014.09.020
14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007
15. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616
16. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. (2009) 4:44–57. doi: 10.1038/nprot.2008.211
17. Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. (2009) 37:1–13. doi: 10.1093/nar/gkn923
18. Bindea G, Galon J, Mlecnik B. CluePedia cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics. (2013) 29:661–3. doi: 10.1093/bioinformatics/btt019
19. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
20. Bader GD, Hogue CWV. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. (2003) 4:2. doi: 10.1186/1471-2105-4-2
21. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and subnetworks from complex interactome. BMC Syst Biol. (2014) 8(Suppl. 4):S11. doi: 10.1186/1752-0509-8-S4-S11
22. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumorinfiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307
23. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337
24. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. (2012) 2:401–4. doi: 10.1158/2159-8290.CD-12-0095
25. Chu VT, Gottardo R, Raftery AE, Bumgarner RE, Yeung KY. MeV+R: using MeV as a graphical user interface for Bioconductor applications in microarray analysis. Genome Biol. (2008) 9:R118. doi: 10.1186/gb-2008-9-7-r118
26. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. (2017) 45:W98–102. doi: 10.1093/nar/gkx247
27. Rooney Michael S, Shukla Sachet A, Wu Catherine J, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. (2015) 160:48–61. doi: 10.1016/j.cell.2014.12.033
28. Ribas A, Wolchok JD. Cancer immunotherapy using checkpoint blockade. Science. (2018) 359:1350–5. doi: 10.1126/science.aar4060
29. Meng X, Huang Z, Teng F, Xing L, Yu J. Predictive biomarkers in PD-1/PD-L1 checkpoint blockade immunotherapy. Cancer Treat Rev. (2015) 41:868–76. doi: 10.1016/j.ctrv.2015.11.001
30. Geoerger B, Kang HJ, Yalon-Oren M, Marshall LV, Vezina C, Pappo A, et al. Pembrolizumab in paediatric patients with advanced melanoma or a PD-L1-positive, advanced, relapsed, or refractory solid tumour or lymphoma (KEYNOTE-051): interim analysis of an open-label, single-arm, phase 1–2 trial. Lancet Oncol. (2019) 21:121–33. doi: 10.1016/S1470-2045(19)30671-0
31. Meng X, Liu Y, Zhang J, Teng F, Xing L, Yu J. PD-1/PD-L1 checkpoint blockades in non-small cell lung cancer: new development and challenges. Cancer Lett. (2017) 405:29–37. doi: 10.1016/j.canlet.2017.06.033
32. Wang X, Guo G, Guan H, Yu Y, Lu J, Yu J. Challenges and potential of PD-1/PD-L1 checkpoint blockade immunotherapy for glioblastoma. J Exp Clin Cancer Res. (2019) 38:87. doi: 10.1186/s13046-019-1085-3
33. Xu F, Jin T, Zhu Y, Dai C. Immune checkpoint therapy in liver cancer. J Exp Clin Cancer Res. (2018) 37:110. doi: 10.1186/s13046-018-0777-4
34. Teng F, Meng X, Kong L, Yu J. Progress and challenges of predictive biomarkers of anti PD-1/PD-L1 immunotherapy: a systematic review. Cancer Lett. (2018) 414:166–73. doi: 10.1016/j.canlet.2017.11.014
35. Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity's roles in cancer suppression and promotion. Science. (2011) 331:1565–70. doi: 10.1126/science.1203486
36. Eil R, Vodnala SK, Clever D, Klebanoff CA, Sukumar M, Pan JH, et al. Ionic immune suppression within the tumour microenvironment limits T cell effector function. Nature. (2016) 537:539–43. doi: 10.1038/nature19364
37. Levental KR, Yu H, Kass L, Lakins JN, Egeblad M, Erler JT, et al. Matrix crosslinking forces tumor progression by enhancing integrin signaling. Cell. (2009) 139:891–906. doi: 10.1016/j.cell.2009.10.027
38. Trimboli AJ, Cantemir-Stone CZ, Li F, Wallace JA, Merchant A, Creasap N, et al. Pten in stromal fibroblasts suppresses mammary epithelial tumors. Nature. (2009) 461:1084–91. doi: 10.1038/nature08486
39. Wang P, Shi Q, Deng WH, Yu J, Zuo T, Mei FC, et al. Relationship between expression of NADPH oxidase 2 and invasion and prognosis of human gastric cancer. World J Gastroenterol. (2015) 21:6271–9. doi: 10.3748/wjg.v21.i20.6271
40. Aoyama T, Hutchins G, Arai T, Sakamaki K, Miyagi Y, Tsuburaya A, et al. Identification of a high- risk subtype of intestinal- type Japanese gastric cancer by quantitative measurement of the luminal tumor proportion. Cancer Med. (2018) 7:4914–23. doi: 10.1002/cam4.1744
41. Chen YL, Fang JH, Lai MD, Shan YS. Depletion of CD4(+)CD25(+) regulatory T cells can promote local immunity to suppress tumor growth in benzo[a]pyrene-induced forestomach carcinoma. World J Gastroenterol. (2008) 14:5797–809. doi: 10.3748/wjg.14.5797
42. Fashoyin-Aje L, Donoghue M, Chen H, He K, Veeraraghavan J, Goldberg KB, et al. FDA approval summary: pembrolizumab for recurrent locally advanced or metastatic gastric or gastroesophageal junction adenocarcinoma expressing PD-L1. Oncologist. (2019) 24:103–9. doi: 10.1634/theoncologist.2018-0221
43. Fuchs CS, Doi T, Jang RW, Muro K, Satoh T, Machado M, et al. Safety and efficacy of pembrolizumab monotherapy in patients with previously treated advanced gastric and gastroesophageal junction cancer phase 2 clinical KEYNOTE-059 trial. JAMA Oncol. (2018) 4:e180013. doi: 10.1001/jamaoncol.2018.0013
44. Shitara K, Ueha S, Shichino S, Aoki H, Ogiwara H, Nakatsura T, et al. First-in-human phase 1 study of IT1208, a defucosylated humanized anti-CD4 depleting antibody, in patients with advanced solid tumors. J Immunother Cancer. (2019) 7:195. doi: 10.1186/s40425-019-0677-y
45. Koemans WJ, Chalabi M, van Sandick JW, van Dieren JM, Kodach LL. Beyond the PD-L1 horizon: in search for a good biomarker to predict success of immunotherapy in gastric and esophageal adenocarcinoma. Cancer Lett. (2019) 442:279–86. doi: 10.1016/j.canlet.2018.11.001
46. Kim ST, Cristescu R, Bass AJ, Kim KM, Odegaard JI, Kim K, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. (2018) 24:1449–58. doi: 10.1038/s41591-018-0101-z
47. Lanier LL, Corliss B, Wu J, Phillips JH. Association of DAP12 with activating CD94/NKG2C NK cell receptors. Immunity. (1998) 8:693–701. doi: 10.1016/S1074-7613(00)80574-9
48. Dietrich J, Cella M, Seiffert M, Buhring HJ, Colonna M. Cutting edge: signal-regulatory protein β 1 Is a DAP12-associated activating receptor expressed in myeloid cells. J Immunol. (2000) 164:9–12. doi: 10.4049/jimmunol.164.1.9
49. Shabo I, Olsson H, Stal O, Svanvik J. Breast cancer expression of DAP12 is associated with skeletal and liver metastases and poor survival. Clin Breast Cancer. (2013) 13:371–7. doi: 10.1016/j.clbc.2013.05.003
50. Cheray M, Bessette B, Lacroix A, Melin C, Jawhari S, Pinet S, et al. KLRC3, a natural killer receptor gene, is a key factor involved in glioblastoma tumourigenesis and aggressiveness. J Cell Mol Med. (2017) 21:244–53. doi: 10.1111/jcmm.12960
51. Chanmee T, Ontong P, Konno K, Itano N. Tumor-associated macrophages as major players in the tumor microenvironment. Cancers. (2014) 6:1670–90. doi: 10.3390/cancers6031670
52. Takamiya R, Ohtsubo K, Takamatsu S, Taniguchi N, Angata T. The interaction between Siglec-15 and tumor-associated sialyl-Tn antigen enhances TGF-β secretion from monocytes/macrophages through the DAP12–Syk pathway. Glycobiology. (2013) 23:178–87. doi: 10.1093/glycob/cws139
53. Yoshida O, Kimura S, Dou L, Matta BM, Yokota S, Ross MA, et al. DAP12 deficiency in liver allografts results in enhanced donor DC migration, augmented effector T cell responses and abrogation of transplant tolerance. Am J Transplant. (2014) 14:1791–805. doi: 10.1111/ajt.12757
54. Kovats S, Gmyrek GB, Akilesh HM, Graham DB, Fuchs A, Yang L, et al. Loss of DAP12 and FcRc drives exaggerated IL-12 production and CD8+ T cell response by CCR2+ Mo-DCs. PLoS ONE. (2013) 8:e76145. doi: 10.1371/journal.pone.0076145
55. Hong Q, Sze CI, Lin SR, Lee MH, He RY, Schultz L, et al. Complement C1q activates tumor suppressor WWOX to induce apoptosis in prostate cancer cells. PLoS ONE. (2009) 4:e5755. doi: 10.1371/journal.pone.0005755
56. Kaur A, Sultan SH, Murugaiah V, Pathan AA, Alhamlan FS, Karteris E, et al. Human c1q induces apoptosis in an Ovarian cancer cell line via tumor necrosis factor pathway. Front Immunol. (2016) 7:599. doi: 10.3389/fimmu.2016.00599
57. Mangogna A, Belmonte B, Agostinis C, Zacchi P, Iacopino DG, Martorana A, et al. Prognostic implications of the complement protein C1q in gliomas. Front Immunol. (2019) 10:2366. doi: 10.3389/fimmu.2019.02366
58. Yamada Y, Arai T, Kojima S, Sugawara S, Kato M, Okato A, et al. Regulation of antitumor miR-144-5p targets oncogenes: direct regulation of syndecan-3 and its clinical significance. Cancer Sci. (2018) 109:2919–36. doi: 10.1111/cas.13722
59. Linnartz-Gerlach B, Bodea LG, Klaus C, Ginolhac A, Halder R, Sinkkonen L, et al. TREM2 triggers microglial density and age-related neuronal loss. Glia. (2019) 67:539–50. doi: 10.1002/glia.23563
60. Takahashi K, Rochford CD, Neumann H. Clearance of apoptotic neurons without inflammation by microglial triggering receptor expressed on myeloid cells-2. J Exp Med. (2005) 201:647–57. doi: 10.1084/jem.20041611
61. Peng Q, Malhotra S, Torchia JA, Kerr WG, Coggeshall KM, Humphrey MB. TREM2- and DAP12-dependent activation of PI3K requires DAP10 and is inhibited by SHIP1. Sci Signal. (2010) 3:ra38. doi: 10.1126/scisignal.2000500
62. Haure-Mirande JV, Wang M, Audrain M, Fanutza T, Kim SH, Heja S, et al. Integrative approach to sporadic Alzheimer's disease: deficiency of TYROBP in cerebral Aβ amyloidosis mouse normalizes clinical phenotype and complement subnetwork molecular pathology without reducing Aβ burden. Mol Psychiatry. (2019) 24:431–46. doi: 10.1038/s41380-018-0255-6
Keywords: gastric cancer, protein-protein interaction, weighted gene co-expression network analysis, biomarker, prognosis
Citation: Jiang J, Ding Y, Wu M, Lyu X, Wang H, Chen Y, Wang H and Teng L (2020) Identification of TYROBP and C1QB as Two Novel Key Genes With Prognostic Value in Gastric Cancer by Network Analysis. Front. Oncol. 10:1765. doi: 10.3389/fonc.2020.01765
Received: 11 February 2020; Accepted: 06 August 2020;
Published: 11 September 2020.
Edited by:
Timothy I. Shaw, St. Jude Children's Research Hospital, United StatesReviewed by:
Jyoti Sharma, Institute of Bioinformatics (IOB), IndiaMaria Paola Paronetto, Foro Italico University of Rome, Italy
Copyright © 2020 Jiang, Ding, Wu, Lyu, Wang, Chen, Wang and Teng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Lisong Teng, lsteng@zju.edu.cn
†These authors have contributed equally to this work