Xiang-Ting Qiu,Yu-Cui Song,Jian Liu,Zhen-Min Wang,Xing Niu,Jing He

Xiang-Ting Qiu,Jian Liu,Zhen-Min Wang,Department of Clinical Laboratory,Linyi Central Hospital,Linyi 276400,Shandong Province,China

Yu-Cui Song,Department of Operating Room,Linyi Central Hospital,Linyi 276400,Shandong Province,China

Xing Niu,Second Clinical College,Shengjing Hospital Affiliated to China Medical University,Shenyang 110004,Liaoning Province,China

Jing He,Department of Pediatric Surgery,Guangzhou Women and Children's Medical Center,Guangzhou Medical University,Guangzhou 510623,Guangdong Province,China

Abstract BACKGROUND Gastric cancer(GC)is the most commonly diagnosed malignancy worldwide.Increasing evidence suggests that it is necessary to further explore genetic and immunological characteristics of GC.AIM To construct an immune-related gene(IRG)signature for accurately predicting the prognosis of patients with GC.METHODS Differentially expressed genes(DEGs)between 375 gastric cancer tissues and 32 normal adjacent tissues were obtained from The Cancer Genome Atlas(TCGA)GDC data portal.Then,differentially expressed IRGs from the ImmPort database were identified for GC.Cox univariate survival analysis was used to screen survival-related IRGs.Differentially expressed survival-related IRGs were considered as hub IRGs.Genetic mutations of hub IRGs were analyzed.Then,hub IRGs were selected to conduct a prognostic signature.Receiver operating characteristic(ROC)curve analysis was used to evaluate the prognostic performance of the signature.The correlation of the signature with clinical features and tumor-infiltrating immune cells was analyzed.RESULTS Among all DEGs,70 hub IRGs were obtained for GC.The deletions and amplifications were the two most common types of genetic mutations of hub IRGs.A prognostic signature was identified,consisting of ten hub IRGs(including S100A12,DEFB126,KAL1,APOH,CGB5,GRP,GLP2R,LGR6,PTGER3,and CTLA4).This prognostic signature could accurately distinguish patients into highand low- risk groups,and overall survival analysis showed that high risk patients had shortened survival time than low risk patients(P <0.0001).The area under curve of the ROC of the signature was 0.761,suggesting that the prognostic signature had a high sensitivity and accuracy.Multivariate regression analysis demonstrated that the prognostic signature could become an independent prognostic predictor for GC after adjustment for other clinical features.Furthermore,we found that the prognostic signature was significantly correlated with macrophage infiltration.CONCLUSION Our study proposed an immune-related prognostic signature for GC,which could help develop treatment strategies for patients with GC in the future.

Key words:Gastric cancer;Immune-related genes;Tumor microenvironment;Immune infiltration;Prognosis;Signature

INTRODUCTION

Gastric cancer(GC)is the fourth most commonly diagnosed malignancy and the second leading cause of cancer-related death worldwide[1].Although the incidence of GC has decreased year by year,the prognosis of GC patients is still not optimistic,especially in China[2].Yet,its pathogenesis remains unclear,therefore,identification of effective biomarkers and therapeutic targets needs to be addressed.

The tumor microenvironment(TME),especially the immune system,plays a pivotal role in the occurrence and development of GC[3].Dysfunction of the immune system assists tumor cells to avoid immune surveillance.Immunotherapy such as programmed death-1(PD-1)blockade has become a promising strategy for GC treatment[4].However,the clinical outcomes of GC patients are still unsatisfactory,and most of novel immunotherapies are still in the early stages of clinical research[5,6].The underlying mechanisms of the immune checkpoint blockade response are complex.Thus,deeper genetic and immunological characterization of GC is required to guide clinicians in selecting and determining the best treatment options.The prognosis of GC is closely related to crosstalk between immune cells and tumor cells[7].Nevertheless,the role of immune-related genes(IRGs)in predicting GC patients’prognosis has not yet been elucidated.

Considering the prognostic potential of IRGs in GC,in this study,we studied immune-related molecular features.We analyzed IRGs using a large amount of transcription data of GC and explored their potential molecular mechanisms.Based on differentially expressed IRGs,we developed an immune-related prognostic model.Our results will help develop treatment strategies for patients with GC.

MATERIALS AND METHODS

GC datasets

Transcriptome RNA-seq data and corresponding clinical information of GC were retrieved from the Genomic Data Commons(GDC)data portal(https://gdc.xenahubs.net),including 375 GC tissues and 32 normal adjacent tissues[8].A list of IRGs were downloaded from the Immunology Database and Analysis Portal(ImmPort)database[9].Furthermore,we derived a list of transcription factors(TFs)from Cistrome database(http://cistrome.org/CistromeCancer/)[10].Mutation data were obtained from the Broad GDAC Firehose.

Differential expression analysis

Differential expression analysis was performed using the edgeR package(http://bioconductor.org/packages/edgeR/)in R[11].The raw data were normalized by Trimmed mean of M values(TMM).The genes with |log fold change(FC)| >1 and false discovery rate(FDR)<0.01 were considered as differentially expressed genes(DEGs).Differentially expressed IRGs and TFs were extracted from these DEGs.

Functional enrichment analysis

The clusterProfiler package was used to annotate DEGs,including gene ontology(GO)and Kyoto Encyclopedia of Genes and Genomes(KEGG)[12].GO terms include biological process,molecular function,and cellular component.P<0.05 was considered to be significantly enriched.

Cox univariate survival analysis

Cox univariate survival analysis of IRGs was performed using the survival package in R.P<0.05 was set as the screening criterion.Survival-related IRGs were identified.Furthermore,differentially expressed survival-related IRGs were considered as hub IRGs.

Protein-protein interaction network

Hub IRGs were analyzed based on the STRING database(https://string-db.org/).Protein-protein interaction(PPI)network was analyzed using Cytoscape.The interactions among proteins could be involved in the progression of diseases.The nodes with high degree could possess potential as hub genes/proteins.

Molecular characteristics

Mutation data of GC samples from Broad GDAC Firehose website were used to analyze genetic alterations of hub IRGs.The results were visualized into waterfall maps using maftools[13].

Construction of an immune-related prognostic signature

Hub IRGs were analyzed by multivariate regression analysis.Then,hub IRGs were selected to conduct a prognostic signature.The risk score was calculated by the expression levels of hub IRGs and Cox regression coefficient.All patients with GC were divided into a high risk group and low risk group according to the median value of risk score.To verify the prognostic potential,the area under the time-dependent receiver operating characteristic(ROC)curve(AUC)was calculated with the survival ROC package[14].The correlation between clinical features and the prognostic signature was evaluated.Immune infiltrate levels of GC were obtained from the TIMER database[15].The tumor-infiltrating immune cells included macrophages,B cells,CD4+T cells,CD8+ T cells,dendritic cells,and neutrophils.The relationships between immune cells and the prognostic signature were calculated.

RESULTS

Identification of DEGs for GC

Four thousand two hundred and fifty-nine DEGs with |logFC| >1 and FDR <0.01 were identified in 375 GC tissues compared to 32 normal adjacent tissues,including 1951 up-regulated and 2356 down-regulated genes,as depicted in volcano plot(Figure 1A).Hierarchical clustering analysis showed that these DEGs can obviously distinguish GC tissue samples from normal adjacent tissue samples(Figure 1B).These DEGs were significantly associated with immune-related biological processes like humoral immune response and acute inflammatory response(Figure 1C).The top ten cellular components and molecular functions enriched by DEGs are depicted in Figure 1D and 1E,respectively.As shown in Figure 1F,the DEGs were mainly enriched in several pathways related with GC,such as neuroactive ligand-receptor interaction,cytokine-cytokine receptor interaction,chemical carcinogenesis and so on.These results suggested that these DEGs might be involved in the pathogenesis of GC.

Identification of differentially expressed IRGs for GC

A list of IRGs were downloaded from the ImmPort database.Differentially expressed IRGs were screened from all DEGs,including 181 up-regulated and 354 downregulated IRGs.The results are visualized into volcano plot(Figure 2A)and heatmap(Figure 2B).In the biological process results,the differentially expressed IRGs were mainly enriched in immune-related processes such as humoral immune response,phagocytosis,and B cell mediated immunity(Figure 2C).As for CC,the genes were in association with immunoglobulin complex,receptor binding,and circulating immunoglobulin complex(Figure 2D).Intriguingly,these IRGs have the molecular immune-related functions like receptor ligand activity,antigen binding,and cytokine activity(Figure 2E).As expected,these differentially expressed IRGs were mainly enriched in immune-related pathways such as cytokine-cytokine receptor interaction,neuroactive ligand-receptor interaction,and viral protein interaction with cytokine and cytokine receptor(Figure 2F).

Identification of differentially expressed TFs for GC

Sixty-seven TFs were differentially expressed between GC tissues and normal tissues,including 37 up-regulated and 30 down-regulated TFs,as shown in volcano plot and heatmap(Figure 3A and B).We further explored their potential functions.The results showed that these TFs are mainly involved in regionalization,pattern specification process,and gland development(Figure 3C).As expected,these TFs could regulate TF complex,nuclear TF complex,and nuclear chromatin(Figure 3D).Similarly,they have the function of TF activity(Figure 3E).According to KEGG pathway enrichment results,these TFs are involved in several cancers(Figure 3F).

Identification of survival-related IRGs for GC

Cox univariate survival analysis of all IRGs was performed using the survival package in R.The results showed that 183 IRGs were significantly related to overall survival of patients with GC(P<0.05).

Identification of hub IRGs for GC

By intersection of differentially expressed IRGs and survival-related IRGs,70 hub IRGs were obtained for GC(Figure 4A).The forest plot of hazard ratios(HRs)suggested that exceptCTLA4,LGR6,andKIR2DS4,other hub IRGs were risk factors for GC(Figure 4B).As shown in PPI network,IL6,F2R,andAGTwere the top three hub genes(Figure 4C).GO enrichment analysis results showed that these hub IRGs were enriched in many biological processes like positive regulation of response to external stimulus,muscle cell proliferation,and peptidyl-tyrosine phosphorylation(Figure 4D).In the cellular component results,receptor complex,vacuolar lumen,and secretory granule lumen were mainly enriched(Figure 4E).As for molecular function,the hub IRGs were significantly associated with receptor ligand activity,growth factor activity,and peptide hormone binding(Figure 4F).As shown in Figure 4G,hub IGRs are mainly involved in GC-related pathways such as neuroactive ligand-receptor interaction,Rap1 signaling pathway,PI3K-Akt signaling pathway,and cytokinecytokine receptor interaction.

Genetic mutations of hub IRGs for GC

We further analyzed the molecular features of hub IRGs for GC.Their genetic alterations were detected.Among 438 samples,323 had genetic alterations.The results showed that deletions and amplifications were the two most common types of genetic mutations(Figure 5).

Figure 1 Identification of differentially expressed genes for gastric cancer.A and B:Volcano plot(A)and heatmap(B)showing differentially expressed genes(DEGs)between gastric cancer tissues and normal tissues.Blue represents down-regulation and red represents up-regulation;C-E:Gene Ontology enrichment results of the top ten DEGs including biological processes(C),cellular component(D),and molecular function(E);F:Kyoto Encyclopedia of Genes and Genomes enrichment results of the top ten DEGs.

Construction of a prognostic signature for GC

Figure 2 Identification of differentially expressed immune-related genes for gastric cancer.A and B:Volcano plot(A)and heatmap(B)showing differentially expressed immune-related gene(IRGs)between gastric cancer tissues and normal tissues.Blue represents down-regulation and red represents upregulation;C-E:Gene Ontology enrichment results of the top ten differentially expressed IRGs including biological processes(C),cellular component(D),and molecular function(E);F:Kyoto Encyclopedia of Genes and Genomes enrichment results of the top ten differentially expressed IRGs.

Figure 3 Identification of differentially expressed transcription factors for gastric cancer.A and B:Volcano plot(A)and heatmap(B)showing differentially expressed transcription factors(TFs)between gastric cancer tissues and normal tissues.Blue represents down-regulation and red represents upregulation;C-E:Gene Ontology enrichment results of the top ten differentially expressed TFs including biological processes(C),cellular component(D),and molecular function(E);F:Kyoto Encyclopedia of Genes and Genomes enrichment results of the top ten differentially expressed TFs.

Figure 4 Identification of hub immune-related genes for gastric cancer.A:The intersection of differentially expressed immune-related gene(IRGs)and survival-related IRGs;B:Forest plot showing prognostic value of hub IRGs.The X-axis represents hazard ratio,and the Y-axis represents the differentially expressed hub IRG;C:Protein-protein interaction network of hub IRGs;D-F:Gene Ontology enrichment results of the top ten differentially expressed transcription factors(TFs)including biological processes(D),cellular component(E),and molecular function(F);G:Kyoto Encyclopedia of Genes and Genomes enrichment results of the top ten differentially expressed TFs.

Multivariate Cox regression analysis was performed based on the hub IRGs using the survival package in R.Ten hub IRGs were selected to construct a prognostic signature for GC.As depicted in the forest plot,each hub IRG in the prognostic signature can accurately predict GC patients’ prognosis(Figure 6A).Among them,S100A12(HR=1.28;95% confidence interval [CI]:1.1-1.48;P=0.001),DEFB126(HR=1.79;95%CI:1.01-3.17;P=0.045),KAL1(HR=1.36;95%CI:1.13-1.63;P=0.001),APOH(HR=1.13;95%CI:1.01-1.28;P=0.034),CGB5(HR=1.2;95%CI:1.04-1.39;P=0.015),GRP(HR=1.33;95%CI:1.06-1.67;P=0.014),andGLP2R(HR=1.45;95%CI:1.17-1.79;P=0.001)were risk factors for GC,whileLGR6(HR=0.86;95%CI:0.78-0.95;P=0.002),PTGER3(HR=0.79;95%CI:0.64-0.99;P=0.037),andCTLA4(HR=0.8;95%CI:0.66-0.97;P=0.021)were protective factors for GC.Risk score was calculated and the patients were divided into a high risk group and low risk group based on the median value of risk score(Figure 6B,C).Heatmap depicts the expression patterns of the ten hub IRGs between the high risk group and low risk group(Figure 6D).Overall survival analysis showed that high risk patients had shortened survival time than t low risk patients(Figure 6E;P<0.0001).The AUC was 0.761,suggesting that the prognostic signature had a high sensitivity and accuracy(Figure 6F).Correlation analysis showed that the prognostic signature was not significantly associated with clinical features including age(Figure 7A),gender(Figure 7B),pathologic T(Figure 7C),pathologic N(Figure 7D),pathologic M(Figure 7E),and pathologic stage(Figure 7F).As shown in multivariate regression analysis results,the prognostic signature could become an independent prognostic predictor after adjustment for other factors including age,gender,pathologic T,pathologic N,pathologic M,and pathologic stage(Table 1).Furthermore,we also found a difference in expression patterns of the ten hub IRGs in the prognostic signature between GC tissues and normal tissues,includingAPOH(Figure 8A;P=0.006),CGB5(Figure 8B;P=0.00097),CTLA4(Figure 8C;P=0.00021),DEFB126(Figure 8D;P=0.096),GLP2R(Figure 8E;P=2e-13),GRP(Figure 8F;P=6.4e-09),KAL1(Figure 8G;P=8.3e-12),LGR6(Figure 8H;P=2.9e-06),PTGER3(Figure 8I;P=3.7e-10),andS100A12(Figure 8J;P=0.0003).

Relationship between immune infiltration and the prognostic signature

To explore whether the hub IRGs in the prognostic signature are involved in the tumor immune microenvironment,the correlation between immune cell infiltration and prognostic signature was analyzed.We found that the prognostic signature was significantly correlated with macrophage infiltration(Figure 9A).However,there was no significant correlation between the prognostic signature and infiltration of other immune cells including B cells(Figure 9B),CD4+ T cells(Figure 9C),CD8+ T cells(Figure 9D),dendritic cells(Figure 9E),and neutrophils(Figure 9F).

Table 1 Univariate and multivariate Cox regression analyses of clinical features and risk score for gastric cancer

DISCUSSION

In this study,we comprehensively analyzed the transcriptome RNA-seq data of GC and identified 4259 DEGs.Among all DEGs,181 up-regulated and 354 down-regulated DEGs were IRGs.These differentially expressed IRGs were mainly enriched in immune-related processes and pathways such as humoral immune response,phagocytosis,B cell mediated immunity,and cytokine-cytokine receptor interaction.Humoral immune response is in association with the progression of GC[16].B cells play a role in regulating the immune response in GC.B cell depletion may be a useful strategy to enhance anti-tumor immune response[17,18].To explore the potential molecular mechanisms of GC,we also analyzed TFs among all DEGs.The results showed that 67 TFs were differentially expressed in GC.Functional enrichment analysis confirmed that these TFs have the function of TF activity.Thus,it is necessary to further explore the potential mechanisms of TFs we identified in the pathogenesis of GC.

One hundred and eighty-three survival-related IRGs were identified for GC by Cox univariate survival analysis.Among them,70 hub IRGs were differentially expressed and associated with overall survival of GC.PPI network analysis indicated thatIL6,F2R,andAGTwere the top three hub genes.For example,a previous study has found that IL-6 could become a target to overcome chemotherapy resistance in GC[19].These hub IGRs are mainly involved in GC-related pathways such as Rap1 signaling pathway[20],PI3K-Akt signaling pathway[21],and cytokine-cytokine receptor interaction[22].Among 438 samples,323 had genetic alterations.Deletions and amplifications frequently occurred for all hub IRGs.Genetic alteration could become a promising target for the therapy of GC[23].In this study,we constructed an immunerelated prognostic signature consisting of ten hub IRGs.Each hub IRG could separately predict the overall survival of patients with GC.Among them,a member of the S100 family,S100A12,as a component of ubiquitinylation complex,could be involved in β-catenin degradation[24].A previous study found that S100A12 could be significantly associated with invasion and metastasis of GC.It has been considered as an independent prognostic factor for GC,which is consistent with our results[25].Furthermore,CGB5 is correlated with a poor prognosis in patients with advanced GC[26,27].It has been reported that LGR6 might be involved in the development of GCviathe PI3K/AKT/mTOR axis[28].However,the functions of other hub IRGs in GC remain unclear.

Our findings suggested that high risk patients had shortened overall survival time than low risk patients based on the median value of risk score.ROC curve analysis confirmed that the prognostic signature had an excellent performance.Our further analysis showed that the prognostic signature could become an independent prognostic factor after adjusting for other prognostic factors including age,gender,pathologic T,pathologic N,pathologic M,and pathologic stage.Moreover,we found that the prognostic signature was positively correlated with the level of macrophage infiltration.M2 macrophages have been shown to be associated with a poor prognosis in a variety of cancers,including GC[29-31].

Figure 5 Mutation landscape of hub immune-related genes for gastric cancer.The horizontal axis represents different samples,and the vertical axis represents immune-related gene.Different colors represent different types of variation.

Figure 6 Construction of a prognostic signature for gastric cancer.A:Forest plot showing the prognostic values of ten hub immune-related gene(IRGs)in the prognostic signature;B:Prognostic index rank;C:Survival status of patients in the two groups;D:Heatmap showing the expression patterns of ten hub IRGs;E:Overall survival analysis results;F:The area under the receiver operating characteristic curve.

Our research has the following limitations.On the one hand,our prognostic markers were based on gene expression profiles.Due to their shortcomings such as high price,long translation cycle,and high requirements for bioinformatics,it is difficult to popularize them in routine clinical applications.However,some alternative methods may be worthy of further exploration,such as screening for optimized signatures from prognostic characteristics through immunohistochemistry.On the other hand,the training cohort of the immune signature model we constructed was from retrospective studies.Therefore,the model requires to be validated by more datasets.

In summary,we have identified differentially expressed IRGs in this study,which may provide a promising perspective for the treatment of GC.We also found that the signature is positively correlated with immune cell infiltration(especially macrophages)and inflammatory responses.The immune gene signature could effectively predict GC patients’ survival,which may be a useful prediction tool to identify patients who will benefit from immunotherapy.

Figure 7 Correlation between clinical features and the prognostic signature for gastric cancer.The clinical features included age(A),gender(B),pathologic T(C),pathologic N(D),pathologic M(E),and pathologic stage(F).

Figure 8 Expression patterns of the ten hub immune-related genes in the prognostic signature between gastric cancer tissues and normal tissues.A:APOH;B:CGB5;C:CTLA4;D:DEFB126;E:GLP2R;F:GRP;G:KAL1;H:LGR6;I:PTGER3;J:S100A12.

Figure 9 Correlation between the prognostic signature and immune cells including macrophages(A),B cells(B),CD4+ T cells(C),CD8+ T cells(D),dendritic cells(E),and neutrophils(F).

ARTICLE HIGHLIGHTS

Research background

Gastric cancer(GC)is the most commonly diagnosed malignancy worldwide.Increasing evidence suggests that it is necessary to further explore genetic and immunological characteristics of GC.

Research motivation

The prognosis of GC is closely related to the crosstalk between immune cells and tumor cells.Nevertheless,the role of immune-related genes in predicting GC patients’prognosis has not yet been elucidated.

Research objectives

In this study,we aimed to construct an immune-related gene signature for accurately predicting the prognosis of patients with GC.

Research methods

Cox univariate survi val analysis was performed to screen survival-related immunerelated genes(IRGs).Differentially expressed survival-related IRGs were considered as hub IRGs.Hub IRGs were selected to conduct a prognostic signature.Receiver operating characteristic(ROC)curve analysis was performed to evaluate its prognostic performance.The correlation of the signature with clinical features and tumorinfiltrating immune cells was analyzed.

Research results

Our study constructed a prognostic signature consisting of ten hub IRGs(includingS100A12,DEFB126,KAL1,APOH,CGB5,GRP,GLP2R,LGR6,PTGER3,andCTLA4),and it could be an independent prognostic predictor for GC.Furthermore,it was significantly associated with immune cell infiltration(especially macrophages).

Research conclusions

We have proposed an immune-related prognostic signature for GC,which may possess prognostic value as a prediction tool for identification of patients who will benefit from immunotherapy.

Research perspectives

The prognostic signature could help develop treatment strategies for patients with GC in the future.