Bioinformatics Analysis Revealing the Correlation between NF-<i>κ</i>B Signaling Pathway and Immune Infiltration in Gastric Cancer

Chen, Qiuxiang; Du, Xiaojing; Ruan, Pengcheng; Ye, Yumei; Zheng, Jin; Hu, Sunkuan

doi:https://doi.org/10.1155/2022/5385456

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2022 | Article ID 5385456 | https://doi.org/10.1155/2022/5385456

Bioinformatics Analysis Revealing the Correlation between NF-κB Signaling Pathway and Immune Infiltration in Gastric Cancer

Qiuxiang Chen,¹Xiaojing Du,²Pengcheng Ruan,³Yumei Ye,¹Jin Zheng,⁴and Sunkuan Hu⁵

Academic Editor: Kamal Shah

Received05 Jan 2022

Revised05 Jul 2022

Accepted13 Jul 2022

Published28 Jul 2022

Abstract

Although the emerging of immunotherapy conferred a new landscape of gastric cancer (GC) treatment, its response rate was of significant individual differences. Insight into GC immune microenviroment may contribute to breaking the dilemma. To this end, the enrichment score of NF-κB signaling pathway was calculated in each GC sample from The Cancer Genome Atlas (TCGA) via ssGSEA algorithm, and its association with immune infiltration was estimated. Based on NF-κB-related genes, a risk score was established and its involvement in immune infiltration, tumor mutational burden (TMB), and N6-methyladenosine (M6A) modification was analyzed in GC. The results showed that NF-κB signaling pathway promoted the infiltration of immune cells in GC. In addition, GC samples were divided into low- and high-risk groups according to a seven-gene (CARD11, CCL21, GADD45B, LBP, RELB, TRAF1, and VCAM1) risk score. Although the high-risk group displayed high immune infiltration and high expression of M6A regulatory genes, it remains in an immunosuppressive microenviroment and whereby suffers a poorer outcome. Of note, most of hub genes were related to immune infiltration and could serve as an independent prognostic biomarker. Conclusively, our study emphasized the crucial role of NF-κB signaling pathway in GC immune microenviroment and provided several candidate genes that may participate in immune infiltration.

1. Introduction

Although the incidence of gastric cancer (GC) decreased in the past decades, it still remains as a major tumor burden all over the world, especially in East Asian regions [1]. Development of the oncotherapy provides more options to GC patients, while patients with advanced GC frequently suffer a tragic outcome. The emerging of immune checkpoint inhibitors (ICIs) is profoundly altering the therapeutic landscape across a spectrum of cancers, including GC [2]. Programmed death-1 (PD1) inhibitor Pembrolizumab is recommended by the 2021 NCCN guideline for the posterior line treatment in GC patients who had , tumor mutational burden- (TMB-) high (TMB-H, ≥10 mutations per megabase), microsatellite instability-high (MSI-H), or deficient mismatch repair (dMMR) [3]. Unfortunately, even under this standard, the response rate to ICIs varies greatly among distinct individuals, due to the fiendishly complicated tumor immune microenvironment. For this reason, insight of the regulatory mechanisms of cancer immunity is of paramount importance to guiding immunotherapy.

NF-κB family consists of five distinct DNA-binding proteins that form various homodimers and heterodimers and thereby drives a series of signaling pathways that could control cell differentiation, proliferation, survival, invasion, angiogenesis during tumorigenesis, and progression [4]. Convincing evidences indicated that NF-κB signaling pathway acts as a driver role during carcinogenesis and progression of GC [5]. It was considered as a potential therapeutic target for GC patients. More recent progress revealed that NF-κB signaling pathway was involved in cancer immune evasion [6, 7]. These studies indicated that NF-κB signaling pathway functions in most of cells in tumor microenvironment, such as tumor-associated macrophages (TAMs), dendritic cells (DCs), myeloid-derived suppressor cells (MDSCs), natural killer (NK) cells, natural killer T (NKT) cells, T cells, and B cells, and regulates the expression of immune checkpoints, such as PD-L1 [4, 6, 7]. However, rare studies investigated the effect of NF-κB signaling pathway on GC immunity, as well as the involving mechanisms.

Interdisciplinary collaboration has been widely applied to the prevention, diagnosis, and therapy in various diseases [8, 9]. For instance, rational mathematical modeling may help in controlling infection diseases, such as the novel corona virus (COVID-19) and hepatitis B [10–12]. The rapid progress of cancer genetics and genomics boosts the establishment of a growing body of public databases, such as The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) [13]. These databases greatly accelerated the application of bioinformatics to the field of cancers. Here, the single-sample gene set enrichment analysis (ssGSEA) was employed to calculate the enrichment score of NF-κB signaling pathway in each GC sample from TCGA, and its correlation with immune infiltration was analyzed in GC. Based on NF-κB-related genes, a seven-gene risk score was further established by using least absolute shrinkage and selection operator (LASSO). Then, the immune infiltration and TMB as well as N6-methyladenosine (M6A) modification were evaluated in different risk score groups. Finally, the association was investigated between the seven genes and prognosis as well as immune infiltration in GC.

2. Materials and Methods

The overall study processes are shown in Figure 1.

2.1. Data Download

GC RNA sequencing (RNA-seq) data (32 normal and 375 tumor) and related clinical data were downloaded from TCGA database, among which 42 samples were excluded for their follow-up time or survival time less than 30 days. GSE62254 () and GSE84437 () with survival data were downloaded from GEO database. The clinical characteristics of GC samples were displayed in Table S1.

2.2. Differentially Expressed Genes and NF-κB-Related Genes

Differentially expressed genes (DEGs) were identified by using “limma” package of R (version 4.1.0) [14]. The cutoff values were set as follows: and . NF-κB related genes () were obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Figure S1) [15].

2.3. Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was used to reveal the underlying pathways [16]. Based on the median enrichment score, GC samples were classified into low (C1) and high (C2) NF-κB signaling pathway groups. The “clusterProfiler” package of R was applied to conduct the biological process (BP) analysis [17]. and were considered as a significant result.

2.4. Gene Set Variation Analysis

Gene set variation analysis (GSVA) was employed to quantify the involvement of BP in each sample using “GSVA” package of the R software [18]. The “limma” package was further employed to calculate the differentially enriched BP between C1 and C2 [14]. The cutoff values were as follows: and .

2.5. Immune Microenvironment Estimation

The immune microenvironment estimation was conducted by ESTIMATE, ssGSEA, MCPcounter, QUANTISEQ, and TIMER algorithms [19–23]. The association between hub genes and immune infiltration was also analyzed by applying TIMER 2.0 database [23]. Tumor immune dysfunction and exclusion (TIDE) algorithm was used to predict the response to ICIs [24].

2.6. Construction and Validation of a Risk Score

According to patients’ survival, LASSO algorithm was employed to establish a risk score in TCGA dataset, which was achieved by the “lars” package of R [25]. Subsequently, samples were divided into low- and high-risk groups. Survival analysis was achieved by the “survival” package. In addition, its prognostic value was validated in GSE62254 and GSE84437 datasets. As for the hub genes, their association with GC survival was estimated by the Kaplan-Meier Plotter online tool [26].

2.7. Somatic Mutation Analysis

To investigate the somatic mutation in GC, the tumor mutational data were downloaded from TCGA and analyzed by the “maftools” package of R [27]. TMB was also calculated in each GC sample from TCGA.

2.8. Relationship between M6A Genes and Risk Score

M6A regulatory genes are comprised of readers (YTHDC1, YTHDC1, YTHDF1, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1, IGF2BP1, IGF2BP2, IGF2BP3, FMR1), writers (METTL3, METTL14, METTL16, WTAP, METL16, RBM15, RBM15B, VIRMA, and ZC3H13) and erasers (ALKBH3, ALKBH5, and FTO) [28]. The association between M6A regulatory genes and risk score was accessed by Pearson correlation analysis and visualized by “ggplot2” package of R.

2.9. Statistical Analysis

The R software (version 4.1.0) or IBM SPSS Statistics 23 was used to conduct statistical analyses. The -test was applicable to normally distributed data and Mann–Whitney test to nonnormally distributed data. Survival time represented the time from diagnosis to the last follow-up or death. Unless otherwise specified, was considered as a statistically significant result.

3. Results

3.1. NF-κB Signaling Pathway Was Significant Association with Immune Infiltration

At first, NF-κB-related genes () were obtained from the KEGG database and ssGSEA algorithm was applied to calculating each sample’s NF-κB signaling pathway score in TCGA dataset. The heatmap showed the expression of these NF-κB-related genes in GC (Figure 2(a)). According to the median value, GC samples were classified into low (C1) and high (C2) NF-κB signaling pathway groups. Compared to C1 type characterized by low enrichment of NF-κB signaling pathway, C2 type had more high-grade samples, more advanced cases, and less adenocarcinoma proportion (Figure 2(b)), indicating that NF-κB signaling pathway was involved in a worse phenotype. The DEGs between C1 and C2 were shown in the volcano plot (Figure 2(c)). Among that, we found that more immune checkpoint-related genes, such as CD274 (PD-L1), CTLA4, and LAG3, enriched in C2 type GC. GSVA was next applied to analyze the differential biology processes between C1 and C2 type GC. The data indicated that the major enriched biology processes were related to cancer immunity (Figure 2(d)). GSEA also demonstrated that these DEGs mainly enriched in immune-related biology processes (Figure 2(e)). These data suggested that the NF-κB signaling pathway was closely concerned in GC’s immune microenviroment, being consistent with previous studies [4, 6, 7].

(a)

(b)

(c)

(d)

(e)

Figure 2

GC samples were clustered by NF-κB signaling pathway score. (a) The expression of NF-κB-related genes in GC samples. (b) GC samples were divided into two cluster (C1 and C2 type) based on NF-κB signaling pathway score (setting median value as cutoff). The bar charts showed the proportions of different grade, T, pathology, and stage in C1 and C2 type. (c) Volcano plot showed the DEGs between C1 and C2 type. The labeled items were immune checkpoint-related genes. (d) Volcano plot showed differential biological processes analyzed by GSVA between C1 and C2 type. (e) GSEA showed that multiple immune-related biological processes were enriched in C2 type. GC: gastric cancer; DEGs: differential expressed genes; GSVA: gene set variation analysis; GSEA: gene set enrichment analysis.

For this reason, we next evaluated the association between NF-κB signaling pathway score and immune infiltration in TCGA dataset. C2 type had a higher immune score, stromal score, and ESTIMATE score, but a lower tumor purity than C1 (Figure 3(a)). The results of QUANTISEQ algorithm showed that C2 had higher B cells (), CD8 T cells (), Tregs (), M1 macrophage (), M2 macrophage (), and neutrophils (), while C1 had higher uncharacterized cells () (Figure 3(b)). Further, TIMER and ssGSEA algorithms found that all estimated immune cells in C2 were significantly higher than that in C1 (Figures 3(c) and 3(d)). These results implied that patients with high level of NF-κB signaling pathway had a high infiltration of immune cells.

(a)

(b)

(c)

(d)

3.2. Establishing a Risk Score Based on NF-κB-Related Genes

Next, a total of 10451 DEGs (9043 upregulated and 1408 downregulated genes) were identified between tumor and normal samples in TGCA dataset (Figure 4(a)). The Venn diagram showed the overlapped genes between these DEGs and NF-κB-related genes, consisting of 31 upregulated and 3 downregulated genes (Figure 4(b)). Based on these genes, LASSO algorithm was used to construct a seven-gene risk score as follows: (Figure 4(c)). Survival analysis revealed that samples with high risk score displayed a poorer overall survival (OS) than those with low risk score (median OS, mOS 779vs. 1407 days, ; Figure 4(d)). The high-risk group also showed a shorter OS than the low-risk group in both validation cohort GSE84437 () and GSE62254 (; Figures 4(e) and 4(f)). The mOS of the high-risk group was 2610 and 1162 days in these two cohorts, respectively, while the mOS of the low-risk group has not reached within the follow-up time. With increasing risk score, the probability of poor prognosis was increased in both discovery and validation queues (Figures 4(g)–4(i)). Followed heatmap showed the expression of these seven hub genes (Figures 4(j)–4(l)). In addition, the high-risk group had a higher level of NF-κB signaling pathway score than the low-risk group (Figure S2). Previous studies have put forward several risk score models in GC [29–32]. These risk score models were further analyzed in TCGA cohort (Table S2), and their efficacy was estimated by receiver operating characteristic curve (ROC) and area under the curve (AUC). The results showed that our seven-gene risk score had the highest AUC value (Figure S3A). Collectively, a NF-κB-related risk score was constructed for predicting GC patients’ prognosis.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 4

Construction of a seven-gene risk score. (a) Volcano plot showed the DEGs between normal and tumor samples in TCGA. (b) Venn plot showed the overlapping genes between DEGs and NF-κB-related genes. (c) The parameters of LASSO algorithm. (d–f). Survival analysis of the risk score in (d) discovery cohort and (e, f) validation cohort. (g–i). Risk score, survival status, survival time in (g) TCGA, (h) GSE84437, and (i) GSE62254 dataset. (j–l). The expression of the seven hub genes in (j) TCGA, (k) GSE84437, and (l) GSE62254 dataset. DEGs: differentially expressed genes; TCGA: The Cancer Genome Atlas; LASSO: least absolute shrinkage and selection operator.

3.3. High-Risk Group Had High Immune Infiltration but Low Response to Immunotherapy

Then, the immune infiltration was estimated between the two risk groups. Obviously, the high-risk group had a higher stromal score (), immune score (), and ESTIMATE score () compared with the low-risk group (Figure 5(a)). Instead, its tumor purity was lower than the low-risk group (). The correlation was estimated between the seven-gene risk score and immune score, stromal score, and ESTIMATE score as well as tumor purity. The results showed that only the seven-gene risk score was significant association with all of the four scores (Figure S3B). QUANTISEQ algorithm indicated that the high-risk group possessed more immune cell infiltration, including B cells (), M2 macrophage (), monocyte (), Tregs (), CD8 T cells (), and nonregulatory CD4 T cells (), but less uncharacterized cells () compared with the low-risk group (Figure 5(b)). MCPcounter algorithm could calculate the absolute abundance of eight immune cells (T cells, NK cells, neutrophils, myeloid dendritic cells, monocytic lineage, cytotoxic lymphocytes, CD T cells, and B lineage) and two stromal cells (fibroblasts and endothelial cells) in heterogeneous tissues from transcriptomic data [21]. Our data showed that all of those cell populations had a more significant infiltration in the high-risk group than in the low-risk group (Figure 5(c)). Further, ssGSEA algorithm reinforced the high level of immune infiltration in the high-risk group (Figure 5(d)). In a word, the high-risk group had a higher infiltration of immune cells than the low-risk group.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5

High-risk group possessed a high immune infiltration and TIDE score. (a–d) The immune infiltration was evaluated by (a) ESTIMATE, (b) QUANTISEQ, (c) MCPcounter, and (d) ssGSEA algorithms between low- and high-risk groups. (e) TIDE score in low- and high-risk groups. (f) The proportion of patients who may response to ICIs in low- and high-risk groups. TIDE: tumor immune dysfunction and exclusion; ssGSEA: single-sample gene set enrichment analysis; immune checkpoint inhibitors: ICIs; , , , and .

Unfortunately, high immune infiltration has not brought a favourable outcome for the high-risk group, implying higher probability of immune evasion. To this end, we evaluated TIDE score—an estimation of two distinct immune evasion mechanisms in tumor, that is dysfunction of tumor infiltration cytotoxic T lymphocytes (CTL) and immunosuppressive factors drove exclusion of CTL, between the two groups [24]. Compared to the low-risk group, the high-risk group had a higher TIDE score (), indicating more probabilities of immune evasion occurred in this group (Figure 5(e)). Of note, the low-risk group had a higher predicted response to ICIs than the high-risk group (47.59% vs. 28.14%, ; Figure 5(f)).

3.4. High-Risk Group Possessed Low Tumor Mutational Burden

Theoretically, the more mutations, the more neoantigens, the higher probabilities for T cell recognition, and the better ICI response [33]. Here, GC mutational data was obtained from TCGA database to evaluate somatic mutations. We found a higher TMB in the low-risk group, which may be the reason for its higher response rate to ICIs (Figures 6(a) and 6(b)). We also investigated the mutation of the seven hub genes in GC and observed that these genes displayed low mutational rate in GC, ranging from 0% to 6% (Figure 6(c)). Then, patients were divided into low- and high-TMB cohorts based on the optimal cutoff value calculated by “survival” package (Figure S3). Results showed that patients with high TMB have a better prognosis than those with low TMB (1686 vs. 801 days, ; Figure 6(d)). Further stratified analysis indicated that the low-risk group has a significant longer OS than the high-risk group in low-TMB (mOS 1407 vs. 640 days; ) cohort, but no significant difference was observed in the high-TMB cohort (; Figures 6(e) and 6(f)), which may be due to the small samples in the high-TMB cohort. These data suggested that low TMB in high-risk score may be an another factor that caused its poor outcome.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6

High-risk group had a low TMB. (a) The waterfall plot showed the top 20 frequent mutations that were occurred in GC. (b) The low-risk group had a higher TMB than the high-risk group. (c) The mutations of the seven hub genes in GC. (d) Survival analysis of TMB in GC. Samples were divided into low- and high-TMB groups according to the optimal cutoff value. (e, f) Stratified analysis showed that the high-risk group tended to a poorer outcome than the low-risk group in low-TMB cohort, but not in high-TMB cohort. TMB: tumor mutational burden; GC: gastric cancer. .

3.5. Association between the Risk Score and M6A

As a reversible epigenetic modification, M6A could affect both messenger RNA and noncoding RNAs in eukaryotes and play crucial roles in diverse cancer pathological processes, including immune evasion [28]. M6A regulatory genes are divided into writers, erasers, and readers, which are often dysregulated in various cancer types [28]. The heatmap showed the expression of common M6A regulatory genes in high- and low-risk groups (Figure 7(a)). The results of Pearson correlation analysis showed that the risk score was positively correlated with the expression of most M6A regulatory genes, such as IGF2BP1, FTO, and ZC3H13 (Figure 7(b)). The expression of the seven hub genes also significantly related to the expression of M6A regulatory genes, and the majority of which was positive correlation. Among that, TRAF2 was most correlated with M6A regulatory genes. Further analysis indicated that the high-risk group had higher expression of METTL16 (), RBM15B (), ZC3H13 (), IGF2BP1 (), PRRC2A (), YTHDF1 (), ALKBH3 (), ALKBH5 (), and FTO () compared to the low-risk group (Figure 7(c)). These results suggested the close implication of NF-κB signaling pathway in M6A regulatory genes.

(a)

(b)

(c)

3.6. Most of Hub Genes Had an Independent Prognostic Significance

Subsequently, we evaluated the potential prognostic value of these seven hub genes in GC through retrieving Kaplan-Meier Plotter online tool. As shown, nearly all the hub genes displayed a prognostic significance for GC patients except VCAM1 (Figure 8(a)). In particular, high expression of CARD11 had a shorter OS than low expression (mOS, 32.6 vs. 93.2 months; ), which was consistent with the predictive significance of other five hub genes, including CCL21 (mOS, 26.8 vs. 36.2 months; ), GADD45B (mOS, 26.3 vs. 35.4 months; ), LBP (mOS, 22.0 vs. 57.6 months; ), RELB (mOS, 23.9 vs. 65.0 months; ), and TRAF1 (mOS, 36.4 vs. 93.2 months; ). In addition, patients with high CARD11 had shorter time of first progression (FP) and postprogression survival (PPS) than those with low CARD11 (Figures 8(b) and 8(c)). Similar results were observed when analyzing CCL21, GADD45B, LBP, RELB, and TRAF1, but not VCAM1. These data indicated that the hub genes other than VCAM1 may be an independent prognostic signature for GC.

(a)

(b)

(c)

3.7. Most of Hub Genes Were Closely Association with Immune Infiltration

Finally, the correlation between these hub genes and immune infiltration was investigated. The results revealed that the hub genes except LBP were positively correlated with the immune score (Figure 9(a)). Further data indicated the significant association of these hub genes with immune cell infiltration (Figure 9(b)). In detail, infiltration of B cells was negatively correlated with LBP, RELB, and TRAF2 but positively correlated with CARD11; infiltration of CD8 T cells, CD4 T cells, and dentritic cells was positively associated with CARD11, CCL21, GADD45B, RELB, and VCAM1; macrophages’ infiltration was positively correlated with CARD11, CCL21, GADD45B, RELB, and VCAM1, while negatively correlated with TRAF2; neutrophils’ infiltration was positively associated with CCL21, GADD45B, RELB, and VCAM1. In summary, these hub genes may play an important role in regulating GC’s immune microenvironment.

(a)

(b)

4. Discussion

In this work, ssGSEA algorithm was applied to calculating a NF-κB signaling pathway score. We found this score was significantly associated with immune infiltration in GC. In addition, a seven-gene risk score was established according to NF-κB-related genes. High risk score group had a higher immune infiltration and M6A level, but a lower TMB, compared to the low risk score group. Further results indicated that the high risk score group tended to an immunosuppressive microenviroment and showed a poor response to ICIs. Finally, most of the hub genes (CARD11, CCL21, GADD45B, LBP, RELB, TRAF1, and VCAM1) had an independent prognostic signature and performed a close connection to immune infiltration in GC.

Deciphering the molecular mechanism of immune evasion in cancer is the lynchpin to achieving the goal of tailored immunotherapy. As abovementioned, NF-κB signaling pathway worked in almost all of the infiltrated cells in tumor microenvironment [34–36]. For instance, NF-κB participates in macrophage polarization and transforms them from a tumor-promoting M2 phenotype to a M1-like cytotoxic phenotype [34]; NF-κB plays an essential role in T cell and B cell activation, as well as development [35, 36]. In GC, IL-1beta activated MDSCs through an IL-1RI/NF-κB pathway, contributing to an immunosuppressive microenvironment, and whereby promoted tumor progress [37]. Our results showed that the high enrichment of NF-κB signaling pathway tends to increase the infiltration of various immune cells. These data indicated that NF-κB signaling pathway played an essential role in GC immune microenvironment. We next constructed a seven-gene risk score based on NF-κB-related genes and divided the samples into low- and high-risk groups with the median as cutoff value. As expected, low- and high-risk groups also displayed distinct immune infiltration. However, high immune infiltration did not bring about a better prognosis in GC. Analyzing the subtype of the infiltrated immune cells, we found that both immunosuppressive (e.g., Tregs, MDSCs, and macrophage M2) and immunostimulative (CD8 T cells, CD4 T cells, and NK cells) cells were significantly enriched in the high-risk group. Of note, the high-risk group possessed a higher TIDE score, indicating more dysfunctional anticancer immune cells and higher possibility of immune evasion in this group [24]. Accordingly, although NF-κB signaling pathway brings more infiltration of immune cells, it may contribute to an immunosuppressive microenvironment in GC, leading to a worse outcome.

Another observation further explained the reason of the poor prognosis in the high-risk group, that is, lower TMB in this group. TMB is considered as a predictive biomarker of response to ICIs, which has been verified retrospectively or prospectively in melanoma, lung cancer, bladder cancer, etc. [38]. Some study revealed that TMB could also be a predictive biomarker for predicting the efficacy of chemotherapy and target therapy [39, 40]. This study agreed that high TMB predicts a better therapeutic response. A systematic pan-cancer analysis estimated the general prognostic impact of TMB in patients with solid tumors based on TCGA database and revealed that the predictive value of TMB varies from different cancer types [41]. Patients with high TMB had a significantly longer OS than those with low TMB () in GC [41], being consistent with our data. Taken together, high TMB in the low-risk group may be another factor for their superior prognosis. In addition, the TIDE algorithm indicated the higher response rate to ICIs in the low-risk group. Although the different TMB may be a major determination factor to this phenomenon, the risk score was also considered as a potential biomarker for predicting response to ICIs in GC.

As a dynamic and reversible methyl-modification, increasing evidences indicated the closely connection between M6A and immune infiltration characterization in various cancers, including GC [42, 43]. M6A modification is mainly regulated by its related genes, which can be classified into readers, writers, and erasers. NF-κB signaling pathway transcriptionally regulated the expression of a broad range of target genes. We found our risk score was significantly associated with several M6A regulatory genes. What is more, the expression of METTL16, RBM15B, ZC3H13, IGF2BP1, PRRC2A, YTHDF1, ALKBH3, ALKBH5, and FTO was significantly upregulated in the high-risk group than in the low-risk group. These genes played a crucial role in tumor microenviroment. For instance, ALKBH5 regulated Tregs and MDSC accumulation via modulating the expression of Mct4/Slc16a3 [44]; FTO facilitated immune invasion and desensitized tumor cells to T cell cytotoxicity [45]; YTHDF1 was correlated with immune cell infiltration but attenuated DCs’ cross-presentation capacity [46, 47]. Therefore, NF-κB signaling pathway may regulate GC’s immune infiltration via affecting M6A modification.

We next focused on the seven NF-κB related genes (CARD11, CCL21, GADD45B, LBP, RELB, TRAF1, and VCAM1). As well known, NF-κB signaling pathway played a crucial role in Helicobacter Pylori-related gastric carcinogenesis and progression. Its related genes, TRAF1, VCAM1, and RELB, also participated in this progress [48–50]. In addition, CCL21 worked as a driver factor via MALAT1/SRSF1/mTOR axis during the progression of GC [51]. Nevertheless, little study revealed the function of CARD11, GADD45B, and LBP in GC. Here, we found that except LBP, these hub genes performed significant correlation with immune infiltration in GC. Previous study also indicated the role of these genes played in immune function [52–54], which has not been reported in GC yet. To this end, study of these genes may deepen our understanding of the mechanism underlying GC’s immune microenviroment.

5. Conclusions

In summary, a clear correlation was revealed between GC immunity and NF-κB signaling pathway, as well as the risk score based on NF-κB-related genes. We also found that NF-κB signaling pathway was significant association with TMB as well M6A level in GC. In addition, most of hub genes performed an independent prognosis value and significant correlation with immune infiltration in GC. These results indicated that NF-κB signaling pathway played a crucial role in GC immunity, and M6A modification may be an important bridge between them. However, only bioinformatics analyses are not enough to clear the mechanisms underlying NF-κB signaling pathway mediated cancer immunity, and further experimental works are necessary. Accordingly, identified hub genes will be the focus in the follow-up study, which may deepen the understanding of cancer immunity and provide a novel strategy for immunotherapy in GC.

Data Availability

All data were obtained from TCGA, GEO, KEGG, Kaplan-Meier online tool, and TIMER 2.0 datasets.

Conflicts of Interest

The authors report no conflicts of interest.

Authors’ Contributions

QXC and SKH provided the overall design. XJD finished the data collection, analysis, and figure layout. QXC, YMY, and PCR contributed to writing manuscript. JZ contributed to the statistical analyses and paper revision. Qiuxiang Chen and Xiaojing Du contributed equally to this work.

Acknowledgments

This work was supported by the Wenzhou Science & Technology Bureau Project (Grant numbers Y20180075 and Y20190070).

Supplementary Materials

Figure S1: the genes enriched in NF-κB signaling pathway. Figure S2: NF-κB signaling pathway score in different risk score groups. Figure S3: comparison of the seven-gene risk score and previous risk score in GC. Figure S4: the cutoff value of TMB in TCGA cohort. Table S1: the clinical characteristics of GC samples in TCGA and GEO. Table S2: the risk score of each model in GC. (Supplementary Materials)

References

F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, and A. Jemal, “Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries,” CA: a Cancer Journal for Clinicians, vol. 68, no. 6, pp. 394–424, 2018.
View at: Publisher Site | Google Scholar
S. S. Joshi and B. D. Badgwell, “Current treatment and recent progress in gastric cancer,” CA: a Cancer Journal for Clinicians, vol. 71, no. 3, pp. 264–279, 2021.
View at: Publisher Site | Google Scholar
National Comprehensive Cancer Network, “NCCN clinical practice guidelines in oncology,” 2021, March 2021, https://www.nccn.org/.
View at: Google Scholar
K. Taniguchi and M. Karin, “NF-κB, inflammation, immunity and cancer: coming of age,” Nature Reviews Immunology, vol. 18, no. 5, pp. 309–324, 2018.
View at: Publisher Site | Google Scholar
S. Chaithongyot, P. Jantaree, O. Sokolova, and M. Naumann, “NF-κB in gastric cancer development and therapy,” Biomedicine, vol. 9, no. 8, p. 870, 2021.
View at: Publisher Site | Google Scholar
B. Eluard, C. Thieblemont, and V. Baud, “NF-κB in the new era of cancer therapy,” Trends in cancer, vol. 6, no. 8, pp. 677–687, 2020.
View at: Publisher Site | Google Scholar
Y. Grinberg-Bleyer, H. Oh, A. Desrichard et al., “NF-κB c-Rel is crucial for the regulatory t cell immune checkpoint in cancer,” Cell, vol. 170, no. 6, pp. 1096–1108.e13, 2017.
View at: Publisher Site | Google Scholar
D. Crosby, S. Bhatia, K. M. Brindle et al., “Early detection of cancer,” Science, vol. 375, no. 6586, article eaay9040, 2022.
View at: Publisher Site | Google Scholar
A. Din, “The stochastic bifurcation analysis and stochastic delayed optimal control for epidemic model with general incidence function,” Chaos, vol. 31, no. 12, p. 123101, 2021.
View at: Publisher Site | Google Scholar
A. Din, Y. Li, T. Khan, and G. Zaman, “Mathematical analysis of spread and control of the novel corona virus (COVID-19) in China,” Chaos, Solitons and Fractals, vol. 141, p. 110286, 2020.
View at: Publisher Site | Google Scholar
A. Din, Y. Li, F. M. Khan, Z. U. Khan, and P. Liu, “On analysis of fractional order mathematical model of hepatitis B using Atangana–Baleanu Caputo (ABC) derivative,” Fractals, vol. 30, no. 1, p. 2240017, 2022.
View at: Publisher Site | Google Scholar
A. Din, Y. Li, A. Yusuf, and A. I. Ali, “Caputo type fractional operator applied to hepatitis b system,” Fractals, vol. 30, no. 1, p. 2240023, 2022.
View at: Publisher Site | Google Scholar
F. Chibon, “Cancer gene expression signatures - the rise and fall?” European Journal of Cancer, vol. 49, no. 8, pp. 2000–2009, 2013.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic acids research, vol. 43, no. 7, p. e47, 2015.
View at: Publisher Site | Google Scholar
M. Kanehisa, Y. Sato, and M. Kawashima, “KEGG mapping tools for uncovering hidden features in biological data,” Protein Science, vol. 31, no. 1, pp. 47–53, 2022.
View at: Publisher Site | Google Scholar
A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences, vol. 102, no. 43, pp. 15545–15550, 2005.
View at: Publisher Site | Google Scholar
T. Wu, E. Hu, S. Xu et al., “clusterProfiler 4.0: a universal enrichment tool for interpreting omics data,” Innovation, vol. 2, no. 3, article 100141, 2021.
View at: Publisher Site | Google Scholar
S. Hänzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, p. 7, 2013.
View at: Publisher Site | Google Scholar
K. Yoshihara, M. Shahmoradgoli, E. Martínez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, no. 1, p. 2612, 2013.
View at: Publisher Site | Google Scholar
D. A. Barbie, P. Tamayo, J. S. Boehm et al., “Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1,” Nature, vol. 462, no. 7269, pp. 108–112, 2009.
View at: Publisher Site | Google Scholar
E. Becht, N. A. Giraldo, L. Lacroix et al., “Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression,” Genome Biology, vol. 17, no. 1, p. 218, 2016.
View at: Publisher Site | Google Scholar
G. Sturm, F. Finotello, F. Petitprez et al., “Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology,” Bioinformatics, vol. 35, no. 14, pp. i436–i445, 2019.
View at: Publisher Site | Google Scholar
T. Li, J. Fu, Z. Zeng et al., “TIMER2.0 for analysis of tumor-infiltrating immune cells,” Nucleic Acids Research, vol. 48, no. W1, pp. W509–W514, 2020.
View at: Publisher Site | Google Scholar
P. Jiang, S. Gu, D. Pan et al., “Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response,” Nature Medicine, vol. 24, no. 10, pp. 1550–1558, 2018.
View at: Publisher Site | Google Scholar
Y. Y. Zhou, L. P. Chen, Y. Zhang et al., “Integrated transcriptomic analysis reveals hub genes involved in diagnosis and prognosis of pancreatic cancer,” Molecular Medicine, vol. 25, no. 1, p. 47, 2019.
View at: Publisher Site | Google Scholar
A. M. Szász, A. Lánczky, Á. Nagy et al., “Cross-validation of survival associated biomarkers in gastric cancer using transcriptomic data of 1,065 patients,” Oncotarget, vol. 7, no. 31, pp. 49322–49333, 2016.
View at: Publisher Site | Google Scholar
A. Mayakonda, D. C. Lin, Y. Assenov, C. Plass, and H. P. Koeffler, “Maftools: efficient and comprehensive analysis of somatic variants in cancer,” Genome Research, vol. 28, no. 11, pp. 1747–1756, 2018.
View at: Publisher Site | Google Scholar
Y. Yang, P. J. Hsu, Y. S. Chen, and Y. G. Yang, “Dynamic transcriptomic m⁶A decoration: writers, erasers, readers and functions in RNA metabolism,” Cell Research, vol. 28, no. 6, pp. 616–624, 2018.
View at: Publisher Site | Google Scholar
M. Yu, Y. Zhang, R. Mao, C. Zhu, R. Zhao, and L. Jin, “A risk model of eight immune-related genes predicting prognostic response to immune therapies for gastric cancer,” Genes (Basel), vol. 13, no. 5, p. 720, 2022.
View at: Publisher Site | Google Scholar
L. Yue, R. Zhang, S. Chen, and G. Duan, “Risk score prediction model of prognosis in gc patients by age and gender combined with m6A modification genes FTO and RBM15,” Frontiers in Cell and Development Biology, vol. 10, p. 710708, 2022.
View at: Publisher Site | Google Scholar
M. Li, W. Cao, B. Huang et al., “Establishment and analysis of an individualized immune-related gene signature for the prognosis of gastric cancer,” Frontiers in surgery, vol. 9, p. 829237, 2022.
View at: Publisher Site | Google Scholar
W. Shao, Z. Yang, Y. Fu et al., “The pyroptosis-related signature predicts prognosis and indicates immune microenvironment infiltration in gastric cancer,” Frontiers in Cell and Development Biology, vol. 9, p. 676485, 2021.
View at: Publisher Site | Google Scholar
D. L. Jardim, A. Goodman, D. de Melo Gagliato, and R. Kurzrock, “The challenges of tumor mutational burden as an immunotherapy biomarker,” Cancer Cell, vol. 39, no. 2, pp. 154–173, 2021.
View at: Publisher Site | Google Scholar
T. Hagemann, T. Lawrence, I. McNeish et al., “"Re-educating" tumor-associated macrophages by targeting NF-kappaB,” The Journal of Experimental Medicine, vol. 205, no. 6, pp. 1261–1268, 2008.
View at: Publisher Site | Google Scholar
S. Gerondakis and U. Siebenlist, “Roles of the NF-kappaB pathway in lymphocyte development and function,” Cold Spring Harbor Perspectives in Biology, vol. 2, no. 5, p. a000182, 2010.
View at: Publisher Site | Google Scholar
S. Gerondakis, T. S. Fulford, N. L. Messina, and R. J. Grumont, “NF-κB control of T cell development,” Nature Immunology, vol. 15, no. 1, pp. 15–25, 2014.
View at: Publisher Site | Google Scholar
S. Tu, G. Bhagat, G. Cui et al., “Overexpression of interleukin-1beta induces gastric inflammation and cancer and mobilizes myeloid-derived suppressor cells in mice,” Cancer Cell, vol. 14, no. 5, pp. 408–419, 2008.
View at: Publisher Site | Google Scholar
A. M. Goodman, S. Kato, L. Bazhenova et al., “Tumor mutational burden as an independent predictor of response to immunotherapy in diverse cancers,” Molecular Cancer Therapeutics, vol. 16, no. 11, pp. 2598–2608, 2017.
View at: Publisher Site | Google Scholar
H. Si, M. Kuziora, K. J. Quinn et al., “A blood-based assay for assessment of tumor mutational burden in first-line metastatic NSCLC treatment: results from the MYSTIC study,” Clinical Cancer Research, vol. 27, no. 6, pp. 1631–1640, 2021.
View at: Publisher Site | Google Scholar
F. Innocenti, F. S. Ou, X. Qu et al., “Mutational analysis of patients with colorectal cancer in CALGB/SWOG 80405 identifies new roles of microsatellite instability and tumor mutational burden for patient outcome,” Journal of Clinical Oncology, vol. 37, no. 14, pp. 1217–1227, 2019.
View at: Publisher Site | Google Scholar
H. X. Wu, Z. X. Wang, Q. Zhao et al., “Tumor mutational and indel burden: a systematic pan-cancer evaluation as prognostic biomarkers,” Annals of translational medicine, vol. 7, no. 22, p. 640, 2019.
View at: Publisher Site | Google Scholar
C. Quan, O. Belaydi, J. Hu et al., “N6-Methyladenosine in cancer immunotherapy: an undervalued therapeutic target,” Frontiers in Immunology, vol. 12, p. 697026, 2021.
View at: Publisher Site | Google Scholar
B. Zhang, Q. Wu, B. Li, D. Wang, L. Wang, and Y. L. Zhou, “m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer,” Molecular Cancer, vol. 19, no. 1, p. 53, 2020.
View at: Publisher Site | Google Scholar
N. Li, Y. Kang, L. Wang et al., “ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment,” Proceedings of the National Academy of Sciences, vol. 117, no. 33, pp. 20159–20170, 2020.
View at: Publisher Site | Google Scholar
R. Su, L. Dong, Y. Li et al., “Targeting FTO suppresses cancer stem cell maintenance and immune evasion,” Cancer Cell, vol. 38, no. 1, pp. 79–96.e11, 2020.
View at: Publisher Site | Google Scholar
J. Hu, D. Qiu, A. Yu et al., “YTHDF1 is a potential pan-cancer biomarker for prognosis and immunotherapy,” Frontiers in Oncology, vol. 11, p. 607224, 2021.
View at: Publisher Site | Google Scholar
D. Han, J. Liu, C. Chen et al., “Anti-tumour immunity controlled through mRNA m⁶A methylation and YTHDF1 in dendritic cells,” Nature, vol. 566, no. 7743, pp. 270–274, 2019.
View at: Publisher Site | Google Scholar
X. K. Wan, S. L. Yuan, H. X. Tao et al., “The upregulation of TRAF1 induced by Helicobacter pylori plays an antiapoptotic effect on the infected cells,” Helicobacter, vol. 21, no. 6, pp. 554–564, 2016.
View at: Publisher Site | Google Scholar
J. Shen, J. Zhai, Q. You et al., “Cancer-associated fibroblasts-derived VCAM1 induced by H. pylori infection facilitates tumor invasion in gastric cancer,” Oncogene, vol. 39, no. 14, pp. 2961–2974, 2020.
View at: Publisher Site | Google Scholar
R. Mejías-Luque, J. Zöller, F. Anderl et al., “Lymphotoxin β receptor signalling executes Helicobacter pylori-driven gastric inflammation in a T4SS-dependent manner,” Gut, vol. 66, no. 8, pp. 1369–1381, 2017.
View at: Publisher Site | Google Scholar
Q. Fu, X. Tan, H. Tang, and J. Liu, “CCL21 activation of the MALAT1/SRSF1/mTOR axis underpins the development of gastric carcinoma,” Journal of Translational Medicine, vol. 19, no. 1, p. 210, 2021.
View at: Publisher Site | Google Scholar
J. M. Salvador, J. D. Brown-Clay, and A. J. Fornace Jr., “Gadd45 in stress signaling, cell cycle control, and apoptosis,” Advances in Experimental Medicine and Biology, vol. 793, pp. 1–19, 2013.
View at: Publisher Site | Google Scholar
M. Lork, J. Staal, and R. Beyaert, “Ubiquitination and phosphorylation of the CARD11-BCL10-MALT1 signalosome in T cells,” Cellular Immunology, vol. 340, p. 103877, 2019.
View at: Publisher Site | Google Scholar
S. Sharma, P. Kadam, and S. Dubinett, “CCL21 programs immune activity in tumor microenvironment,” Advances in Experimental Medicine and Biology, vol. 1231, pp. 67–78, 2020.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Qiuxiang Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

177

Downloads

310

Citations