Molecular Analysis of Prognosis and Immune Pathways of Pancreatic Cancer Based on TNF Family Members

Background Tumor necrosis factor (TNF) family members play a vital role in anticancer therapy. This study aimed to screen the critical markers for the prognostic analysis of pancreatic adenocarcinoma (PAAD) by analyzing the clustering patterns of TNF family members in PAAD. Methods In this study, the NMF clustering method was adopted to cluster samples from The Cancer Genome Atlas (TCGA) to acquire the clustering pattern of the TNF family in PAAD. Differential gene analysis was performed according to TNF family gene clusters. The support vector machine (SVM) method was conducted for further gene screening, and the risk score model of the screened genes was constructed by Lasso. The single sample gene set enrichment analysis (ssGSEA) method was adopted for immunoenrichment analysis and tumor immune cycle analysis. Genes associated with risk scores were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Results We clustered PAAD into two groups based on TNF family genes. Nineteen TNF family genes were significantly associated with the clinical characteristics of PAAD patients. The risk score formula was composed of RHOD, UBE2C, KLHDC7b, MSLN, ADAM8, NME3, GNG2, and MCOLN3. GSE57495 and GSE62452 datasets verified that patients in the high-risk group had a worse prognosis than those in the low-risk group. The risk score-related genes analyzed by GO and KEGG were mainly involved in the modulation of chemical synaptic transmission and synaptic vesicle cycle pathway. There were significant differences in the expression of 15 immune cells between the high-risk group and the low-risk group. The risk score was positively correlated with HCK, interferon, MHC-I, and STAT1. The expression of genes relevant to chemokine, immunostimulator, MHC, and receptor was strongly associated with the risk score. Conclusion The risk score model based on the TNF family can predict the prognosis and immune status of PAAD patients. Further research is needed to verify the clinical prognostic value of risk scores.


Introduction
Pancreatic adenocarcinoma (PAAD) has a 5-year survival rate of less than 10%, which is one of the most aggressive malignant tumors with the worst prognosis [1,2]. About 90% of PAADs are ductal adenocarcinomas originating from the glandular epithelium. In recent years, the morbidity and mortality of the disease have increased significantly. Studies have predicted that PAAD may become the second leading cause of cancer-related death after 2030 [3]. Low immunophenotype and low tumor sensitivity to cytotoxic drugs are the primary causes for the low survival rate in PAAD patients [4]. e cytotoxic chemotherapy for PAAD is primarily limited to the absence of molecular markers to predict the efficacy of chemotherapy [5]. e lack of a specific tumor prognostic model is a serious challenge in the treatment of PAAD. erefore, more efforts are urgently needed to understand the specific prognostic markers and immune function of PAAD.
TNF and TNF receptor (TNFR) superfamily (TNFSF/ TNFRSF) consist of nineteen ligands and twenty-nine receptors [6]. TNF family members are expressed naturally by the immune system and kill tumor activity [7]. In addition, tumor necrosis factor-(TNF-) related ligands have high expectations in anticancer therapy due to their induction of apoptosis [8]. For example, tumor necrosis factor-related apoptosis-inducing ligand (TRAIL), a member of the TNF family, has been shown to selectively induce apoptosis in cancer cells by binding or trimerizing its functional receptors [9]. In PAAD tissues, tumor-infiltrated immature M0 macrophages exhibit antitumor activity by secreting TNF-α [10]. Although a large body of evidence has shown that TNF family members play an essential role in many cancers, including PAAD, the role of TNF family members in PAAD is still not systematically understood.
In recent years, bioinformatics analysis of biodata information obtained from the public databases has played a very positive role in better understanding and treatment of several diseases, including cancer. Using algorithms such as ssGSEA, the expression level of markers can reflect the infiltration of specific cell types in tumor tissue. rough the complete follow-up data of multiple cohorts, the relevance between the relative infiltration level of specific cell types and the survival rate of patients can be determined [11]. For instance, based on the epigenetic properties of immunomodulatory cytokine genes, methylation of these genes has been found to be related to overall survival (OS), disease-specific survival, and disease progression in PAAD patients [12]. Yao et al. discovered differential splicing of AS events between PAAD and normal tissues and successfully constructed a prognostic model to predict the prognosis of PAAD patients using survival-related splicing factors [13]. More recently, Zhang et al. constructed a prognostic model based on the TNF family to predict the prognosis and immune status of lung adenocarcinoma (LUAD) patients [14]. In addition, it has been found that ubiquitin-specific protease 4 (USP4) plays a tumor-promoting role in PAAD and can be used as a prognostic indicator and therapeutic target in patients with PAAD resection [15]. However, no details of expression of TNF family members in PAAD and their clinical significance have been reported.
is study was a systematic study of expression patterns of TNF family members and their clinical significance in PAAD. We aimed to establish a prognostic model for PAAD based on the TNF family by in-depth analysis of relevant data from the TCGA and gene expression integration (GEO) databases. We hope that the prognostic risk score of this study will contribute to the prognosis of PAAD and the formulation of phase immunotherapy strategies.

Datasets and Preprocessing.
e TCGA pancreatic adenocarcinoma dataset (TCGA-PAAD) was downloaded from the TCGA website, and 178 samples were included in the analysis. GSE57495 and GSE62452 datasets were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). Affymetrix generated raw data from the microarray dataset. e RMA algorithm in the Affy package was then applied to perform quantile normalization and background correction for Affymetrix raw data. GSE57495 included 63 samples, and GSE62452 included 66 samples.

Clustering
Based on the TNF Family. TNF family was obtained as per the previous paper [14]. TCGA-PAAD was clustered by the NMF clustering method, and the clustering pattern based on the TNF family was obtained.

Establishment of the Risk Score Model.
According to the TNF family gene cluster, the differential gene analysis was performed for the two clusters. e differential gene of standard was defined as |logFC| > log2(1.5) with P < 0.05. e independent prognostic significance of TNF family members was assessed by univariate Cox regression analysis (P < 0.05). HR > 1 was considered a prognostic risk gene, while HR < 1 was considered a protective prognostic gene. Univariate analysis was carried out for differential genes, and then SVM was used for further screening. e selected genes were modeled using Lasso, and the risk score was the sum of gene expression value * Lasso coefficient.

Immunoinfiltration Analysis.
e ssGSEA method was used for immunoinfiltration analysis. e expression levels of twenty-eight types of cells were mainly analyzed, including immature dendritic cell, immature B cell, activated B cell, activated CD4 T cell, activated CD8 T cell, macrophage, mast cell, MDSC, memory B cell, monocyte, activated dendritic cell, CD56bright natural killer cell, CD56dim natural killer cell, central memory CD4 T cell, central memory CD8 T cell, effector memory CD4 T cell, effector memory CD8 T cell, natural killer cell, natural killer T cell, neutrophil, plasmacytoid dendritic cell, regulatory T cell, follicular T-helper cell, type 1 T-helper cell, type 17 T-helper cell, type 2 T-helper cell, eosinophil, and gamma delta T cell [16]. Tumor immune-cycle analysis mainly analyzed seven steps of immune activity [17].

Pathway Analysis.
Correlation analysis was performed on risk score and all genes, and the correlation standard was defined as |cor| > 0.6. Related genes were analyzed for functional enrichment, mainly by GO and KEGG analysis.

Statistical Analysis.
e R package SurvMiner was used to draw all survival curves. e normality of variables was checked by the Shapiro-Wilk normality test. e unpaired Student's t-test was used to compare the differences between the two groups for variables that conform to the normal distribution. e nonnormally distributed variables were analyzed by the Wilcoxon test. e Kaplan-Meier method was used to generate and visualize subgroup survival curves. e logarithmic rank test was used to determine the statistical significance of the differences in each dataset. All heat maps were generated based on PHEATMAP. All statistical analyses were performed in R (https://www.r-project.org/, version 3.5.1). All the tests were two-sided, and P values < 0.05 were considered statistically significant.

Functional Analysis of the Prognostic Model.
To understand the function of the obtained risk score-related genes, we first conducted correlation analysis of risk scorerelated genes and clinical characteristics, as shown in Figure 4(a). To further understand the potential functions of these genes, GO and KEGG enrichment analyses were performed, respectively. As shown in Figure 4(b), GO analysis revealed that these genes mainly participate in the modulation of chemical synaptic transmission, regulation of trans-synaptic signaling, neurotransmitter secretion, signal release from synapse, regulation of exocytosis, calcium ionregulated exocytosis, regulation of membrane potential, vesicle-mediated transport in synapse, signal release, regulation of calcium ion-dependent exocytosis, regulation of synaptic plasticity, synaptic vesicle transport, establishment of synaptic vesicle localization, synaptic vesicle exocytosis, and synaptic vesicle cycle. KEGG analysis presented that these genes were primarily connected with synaptic vesicle cycle, insulin secretion, and dopaminergic synapse. ese results manifest that the functions of these genes are mainly embodied in the regulation of information transmission between cells and the transfer of nanoparticles. In addition, we conducted GSEA to comprehensively define the features of risk score. e hallmark gene set enrichment analysis showed that risk score-related genes were enriched in MYC targets V2 and TGF-β signaling ( Figure S5A). GO enrichment analysis found that the genes were enriched in pore complex assembly and cysteine-type endopeptidase activity involved in the apoptotic process ( Figure S5B). KEGG enrichment analysis revealed that the genes were concentrated in pathways involving pentose phosphate pathway ( Figure S5C).

Immune Infiltration and Inflammation Analysis.
It is well known that immune cell infiltration is closely related to inflammation. To this end, the ssGSEA method was used for immunoenrichment analysis of 28 kinds of immune cells. As shown in Figure 5(a), the expression of 15 kinds of immune cells showed a notable difference between the high-risk group and low-risk group, including activated B cell, activated CD8 T cell, CD56dim natural killer cell, effector memory CD4 T cell, effector memory CD8 T cell, eosinophil, immature B cell, immature dendritic cell, macrophage, mast cell, MDSC, monocyte, plasmacytoid dendritic cell, T follicular helper cell, and type 1 T-helper cell. en, the immune-cycle activity scores of the two groups were statistically analyzed. We noted that Step2 (cancer antigen presentation), Step4 (including CD4 T cell, dendritic cell, macrophage, T cell, TH17 cell, and Treg cell recruitment), and Step6 (recognition of cancer cell by T cell) displayed statistical differences in anticancer immunity between the two groups ( Figure 5(b)). As the risk score changed, the expression of genes in different inflammatory marker gene sets also changed accordingly, as shown in Figure 5(c). Genome set variation analysis (GSVA) was used to analyze the results of these 7 gene classes: HCK, IgG, interferon, LCK, MHC-I, MHC-II, and STAT1. After analyzing the correlation between the risk score and inflammatory indicators, we noticed that the risk score was significantly positively correlated with HCK, interferon, MHC-I, and STAT1 (P < 0.05) ( Figure S2).

Prediction of the Risk Score for Immunotherapy Response.
We used the TIDE algorithm to verify the risk score of the anti-PD-1 immunotherapy cohorts IMvigor 210 and GSE78220. In the GSE78220 cohort, the risk score in the complete/partial response (CR/PR) group was low compared to the stable/progressive disease (SD/PD) group (P � 0.04, Figure S6A). e CR/PR group had a higher percentage of scores than in the SD/PD group in the IMvigor 210 cohort (Figures S6B and S6C). Besides, the GSE79668 cohort confirmed that the prognosis of the high-risk group was worse than that of the low-risk group ( Figure S6D).
We also analyzed the correlation between risk scores and classical immune checkpoint molecules based on antigen present, cell adhesion, coinhibitor, costimulator, ligand, receptor, and other classifications. We found that the risk score was positively correlated with MICA, ICAM1, CD276, CD80, and TNFSF9 (P < 0.05, Figures S7A-E), and negatively correlated with ADORA2A and AEG1 (P < 0.05, Figures S7F-G).

Discussion
rough bioinformatics analysis, we noted that the TNF family genes were significantly associated with the clinical characteristics of PAAD patients. Our prognostic model based on TNF family genes has been proved clinically adaptable in predicting OS in PAAD patients. e risk score is an independent prognostic risk factor. In the PAAD prognostic model, the expression of immune cells, immunecycle activity, and inflammatory markers was closely correlated with risk assessment. ere was also obvious relativity between risk scores and immune checkpoints. e tumor microenvironment of PAAD is highly immunosuppressive [18]. e complex tumor microenvironment has become one of the challenges that impedes PAAD treatment and leads to immune escape of pancreatic malignant cells [19]. TNF is not only a pleiotropic cytokine that triggers NF-κB activation or RIPK1 kinase-dependent cell death but also a major mediator in inflammation [20]. As a type II transmembrane protein, TNF-α binds to tumor necrosis factor receptor 1 (TNF-R1) and TNF-R2, which subsequently activates downstream signaling pathways [7,21]. Interestingly, TNF plays a "double-edged sword" role in cancer, largely depending on the role of TNF-R1 and TNF-R2 [22,23]. TNF acts as a cancer suppressor by binding to TNF-R1. TNF-R2, on the other hand, can transform the tumor-suppressive TNF into the tumor promoter [24]. TNFα and TGF-β, members of the TNF family, also play an important role in regulating TME [25,26]. Low levels of TNF-α could increase tumor growth by inducing recruitment of endothelial phenotypes of monocytes to the tumor site [27]. Some researchers have proposed that local enhancement of endogenous TNF-α activity can accelerate the death of tumor cells without the associated systemic toxicity [28]. TNF/TNFR superfamily proteins are the major regulatory factors of T cells, among which Fas, TNF-R1, and TRAILR play an important role in promoting apoptosis and inhibiting T-cell activity [29]. Besides, TNFSF10 polymorphism has been identified as a possible prognostic factor for survival in patients undergoing surgery for invasive breast cancer [30]. ese encouraging studies hint at the potential of TNF family members in their efforts to diagnose and predict cancer. In this study, we clustered PAAD into two classes according to TNF family genes, acquired nineteen TNF family genes that were significantly correlated with the clinical characteristics (M, N, T, stage, grade, gender, and age) of PAAD patients, and further screened and established a risk score model. is model has been shown to have certain clinical prognostic value in PAAD.
Stimulation or inhibition of TNF superfamily signaling pathways may influence tumor progression [31]. We compared our model with the model established by Rong et al. [32]. In terms of sample size, our study included more samples (n � 178). Moreover, based on the time-dependent ROC results, the AUC of our model at 1 year, 3 years, and 5 years are 0.714, 0.794, and 0.844, relatively (Figure 3(c)). e 1-year, 3-year, and 5-year AUCs in the model established previously were 0.707, 0.75, and 0.795, respectively ( Figure  S4). is means that our model has better adaptability. We also noticed an interesting finding that the sensitivity of risk scores to predict 1-year, 3-year, and 5-year survival rates gradually increased. e risk factors included in our model, such as UBE2C and ADAM8, and their high expression levels are associated with poor clinical outcomes [33,34]. erefore, we speculate that the increased sensitivity of the risk score to survival may be due to the stronger expression of related genes with the development of PAAD.
Immune cells account for nearly 50% of the components of pancreatic ductal adenocarcinoma cells, but only a few are antitumor effector cells [35]. A novel study pointed out that monocytes may be an effective predictor of response to treatment in patients with glioma [36]. Our research also noticed that monocytes were significantly different between the high-risk group and low-risk group. Immunoinfiltration of macrophages has been used as a prognostic factor to assess the immune microenvironment of pancreatic ductal carcinoma [37]. is study found a significant correlation between the TNF family-related risk score model and the immune infiltration and inflammatory indicators. Immune checkpoint blockade (ICB) therapy is one of the most promising immunotherapies, especially in inhibiting metastasis. However, due to the immunosuppressive tumor microenvironment and extensive fibrotic matrix, immunotherapy is still greatly hindered in the treatment of pancreatic cancer [38]. For example, ICBs with programmed cell death protein-1 (PD-1)/programmed cell death ligand 1 (PD-L1) antibodies showed a sustained response rates in immunogenic tumors [39]. Patients with a large number of tumor-infiltrating lymphocytes in pancreatic ductal adenocarcinoma after ICB apparently have a better prognosis, while patients with mismatched repair defects have a better outcome, suggesting the possibility of a comprehensive immune enhancement that reverses the tumor microenvironment [40]. rough the analysis of different categories of immune checkpoints, we found that immune checkpoints changed as risk scores changed from low to high, and there was a strong correlation between the two. Zhang et al. proposed that the exploration of more valuable PD-L1 and CTLA-4 modulators to improve the efficacy of immunotherapy is currently an effective strategy to promote personalized cancer therapy [41].
In conclusion, the risk score model based on TNF family has good clinical value and adaptability in the prognosis of PAAD. Patients with high-risk PAAD tend to have a poorer prognosis, particularly with respect to immune infiltration and inflammation. In addition, we noted a strong correlation between risk scores and immune checkpoints. Our study may provide novel guidance for the diagnosis, prognosis, and treatment of PAAD.

Conflicts of Interest
e authors declare no conflicts of interest.