Construction and Validation of an Immune-Based Prognostic Model for Pancreatic Adenocarcinoma Based on Public Databases

Background: Pancreatic adenocarcinoma (PAAD) is a highly lethal and aggressive tumor with poor prognoses. The predictive capability of immune-related genes (IRGs) in PAAD has yet to be explored. We aimed to explore prognostic-related immune genes and develop a prediction model for indicating prognosis in PAAD.

Methods: The messenger (m)RNA expression profiles acquired from public databases were comprehensively integrated and differentially expressed genes were identified. Univariate analysis was utilized to identify IRGs that related to overall survival. Whereafter, a multigene signature in the Cancer Genome Atlas cohort was established based on the least absolute shrinkage and selection operator (LASSO) Cox regression analysis. Moreover, a transcription factors regulatory network was constructed to reveal potential molecular processes in PAAD. PAAD datasets downloaded from the Gene Expression Omnibus database were applied for the validations. Finally, correlation analysis between the prognostic model and immunocyte infiltration was investigated.

Results: Totally, 446 differentially expressed immune-related genes were screened in PAAD tissues and normal tissues, of which 43 IRGs were significantly related to the overall survival of PAAD patients. An immune-based prognostic model was developed, which contained eight IRGs. Univariate and multivariate Cox regression revealed that the risk score model was an independent prognostic indicator in PAAD (HR > 1, P < 0.001). Besides, the sensitivity of the model was evaluated by the receiver operating characteristic curve analysis. Finally, immunocyte infiltration analysis revealed that the eight-gene signature possibly played a pivotal role in the status of the PAAD immune microenvironment.

Conclusion: A novel prognostic model based on immune genes may serve to characterize the immune microenvironment and provide a basis for PAAD immunotherapy.

Introduction

Pancreatic cancer (PC) is one of the deadliest and most aggressive malignant neoplasms worldwide (Ilic and Ilic, 2016). In the next decade, PC is estimated to be the second leading cause of death among malignant cancer-related diseases (Rahib et al., 2014; Ferlay et al., 2016). Pancreatic adenocarcinoma (PAAD) occurs in approximately 85% of all PC cases and is associated with a less favorable prognosis (Higuera et al., 2016).

Pancreatic cancer treatment comprises surgery, chemotherapy, radiotherapy, neoadjuvant therapy, targeted molecular therapy, and immunotherapy. Nevertheless, the therapeutic effect of these strategies for PAAD is limited. Therefore, accurate prediction of the prognosis can determine if the patient will benefit from more radical treatment, thereby providing the patient with “individualized” systemic treatment to improve the prognosis.

Pancreatic adenocarcinoma is characterized by the high complexity of stomal tissue, which includes immune cells, various growth factors, the extracellular matrix, and fibroblasts. The tumor microenvironment (TME) accounts for about 15–85% of the whole tumor component in PAAD (Erkan et al., 2012; Liang et al., 2017). The complex and heterogeneous TME induced by interactions between pancreatic epithelial/cancer cells and stromal cells is responsible for PC progression and has been implicated in resistance to chemotherapy and immunotherapy (Markowitz et al., 2015; Incio et al., 2016; Ren et al., 2018). Besides, components of the PAAD microenvironment that contribute to immunosuppression correlate with a poor prognosis of patients (Tang et al., 2014; Whatcott et al., 2015; Wang et al., 2017). With deepening of the understanding of the microenvironment of PC, TME-based clinical and translational therapies could be a breakthrough hotspot in PC treatment in the future.

With the remarkable progress of bioinformatics analysis, in many studies, the mining of public databases has been used increasingly to predict cancer prognosis. Among them, immune-related genes (IRGs) have shown an increasingly prominent role in cancer development and immunotherapy (Ge et al., 2019; Huang et al., 2020; Kong et al., 2020; Yang et al., 2020). Predictive biomarkers related to the tumor immune microenvironment are expected to identify additional target molecules and to enhance immunotherapy efficacy (Taube et al., 2018; Bianco et al., 2019; Jiang et al., 2019; Liu et al., 2020; Zhao B. et al., 2020; Zwing et al., 2020). Currently, PC still lacks prognostic biomarkers related to the tumor immune microenvironment. Therefore, it is necessary to explore important biomarkers in PAAD to guide appropriate treatment options to improve the therapeutic efficacy of patients.

In our research, we investigated the messenger (m)RNA expression and corresponding clinical information of PAAD patients from public databases. Next, we constructed an IRGs-based prognostic model in The Cancer Genome Atlas (TCGA) cohort and validated it in the Gene Expression Omnibus (GEO) dataset. The regulatory network structured by differentially expressed transcription factors (DETFs) and prognosis-related IRGs may provide a theoretical basis to reveal the potential mechanisms at the molecular level. Finally, analyses of prognostic “gene signatures” and infiltration of immune cells may provide new ideas for the role of IRGs in predicting PAAD prognosis.

Materials and Methods

Data Acquisition

The transcriptome sequencing data and corresponding clinical data of 176 PAAD patients were extracted from TCGA (172 PAAD specimens and four normal tissue specimens). The RNA-sequencing data of normal pancreatic tissue were acquired from the Genotype-Tissue Expression (GTEx) Project¹ as well (Carithers et al., 2015). It contains the RNA-expression profile of 167 normal pancreatic tissues. Meanwhile, RNA sequencing fragments per kilobase of exon model per million reads mapped (FPKM) data were also obtained for further analyses. For validation cohort, gene expression matrix files and clinical data of 125 patients with PAAD in the GSE71729 dataset were downloaded from the GEO². Match the gene symbols corresponding to the probes according to the annotation file provided by the manufacturer. If a single gene matches multiple probes, the median ranking value accounts for the expression value. We normalized gene expression value using the robust multiarray average (RMA) algorithm, and the normalized data were log2-transformed for further analyses. Publishing guidelines provided by the GEO database were observed, Therefore, there was no requirement for additional ethical approval. Furthermore, a list of IRGs was acquired from the Immunology Database and Analysis Portal (ImmPort) database that shares resources for immunology-related research³ (Bhattacharya et al., 2014). Then, a list of transcription factors (TFs) was obtained from the Cistrome Project⁴, including 318 TFs (Mei et al., 2017).

Analyses of Differentially Expressed Genes in PAAD

The “limma” R package⁵ (Ritchie et al., 2015) was used for analyses of differential expression. Differential gene expression was defined with adjusted-P < 0.01 and | log2 fold change| > 2 as the cutoff criteria. Then, we extracted differentially expressed immune-related genes (DEIRGs) and DETFs from all DEGs based on the lists obtained from ImmPort and Cistrome Cancer databases.

Analyses of DEIRGs in PAAD Using the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Databases

The functions and pathway enrichment of candidate DEIRGs were analyzed using Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8⁶ (Dennis et al., 2003). To explore the underlying biological functions of DEIRGs, the GO and KEGG databases were searched using the R packages “GOplot⁷” (Walter et al., 2015) and “clusterProfiler⁸” (Kanehisa et al., 2017), respectively. Moreover, the cutoff value for pathway screening and significant functionality was placed at P < 0.05. To explain the correlation between enriched pathways and prognostic IRGs, an interaction network was constructed for visual representations.

Transcription Factors-Mediated Prognosis-Related IRGs Modulation Network

A short duration of follow-up usually limits the accuracy of survival analyses. Hence, we selected patients whose duration of follow-up was ≥60 days. To investigate the prognosis-related DEIRGs in PAAD patients, the “survival” R package⁹ was applied to implement the univariate Cox regression analysis (P < 0.01). To explore the interactions between DETFs and prognosis-related DEIRGs, the correlation test function was employed (set thresholds: P < 0.001 and correlation coefficient > 0.4).

Construction of the IRGs-Based Prognostic Model for PAAD

An IRG prognostic model was developed based on the LASSO Cox regression analysis. To minimize the risk of overfitting, the Lasso method used 10-fold cross validation based on the “glmnet” package¹⁰ in R (R Project for Statistical Computing, Vienna, Austria) (Tibshirani, 1997; Simon et al., 2011). Then, we used β coefficients of the LASSO Cox regression analysis to establish the DEIRGs-based prognostic model for PAAD. We used it to establish a formula to predict the risk score of each patient. The receiver operating characteristic (ROC) curve was used to judge the discrimination ability of various statistical methods on the basis of the binary gold standard (Hanley and McNeil, 1982). The ROC curve was created by the “survival ROC” R package¹¹ to evaluate the sensitivity of the model. Finally, principal component analysis (PCA) was done based on the “prcomp” function from the “stats” R library.

Correlation Between the Immune-Related Signature and Clinical Features in a Prognostic Model of PAAD

The relevance between clinical characteristics (age, gender, histology grade, tumor stage, T staging, N staging, M staging, residual tumor, and outcomes) and expression of eight prognosis signatures in the prognostic model was analyzed using the “beeswarm” R package.

Further Verification of a Prognostic IRG Signature

To verify the prognostic value of the immune-related signature risk score model, we used the GSE71729 dataset as the validation cohort. Samples in the GSE71729 cohort were then divided into high-risk and low-risk groups based on the optimal cut-off point. Kaplan–Meier and ROC curve analysis of the eight-gene signature were performed as mentioned above. In addition, the Human Protein Atlas¹² (Pontén et al., 2011) was used to extract the protein expression of prognostic-related immune genes in tumor samples and normal samples.

Analysis of Immunocyte Infiltration

The Tumor Immune Estimation Resource (TIMER) was employed to analyze and visualize the abundance of tumor-infiltrating immunocytes¹³ (Li et al., 2017). It detailed the abundance of six subsets of tumor-infiltrating immunocytes: B cells, CD8⁺T cells, CD4⁺T cells, macrophages, neutrophils, and dendritic cells (DCs). The online “Immune Estimation” file was retrieved, and the potential correlation between the prediction model and tumor-infiltrating immunocytes was conducted in R.

Statistical Analysis

Statistical analysis was undertaken with R v3.6.3. Unless specified otherwise, P < 0.05 was considered significant.

Results

The flowchart of our study was displayed in Figure 1. The clinical features of the 185 PAAD patients enrolled in the TCGA–PAAD cohort were presented in Table 1.

FIGURE 1

Figure 1. Study flowchart.

TABLE 1

Table 1. Clinical features of patients with pancreatic adenocarcinoma (PAAD).

Identification of DEIRGs and TFs in PAAD

A total of 343 tissues were analyzed [172 PAAD tissues and 171 normal tissues (167 from the GTEx database)]. Compared with normal tissue specimens, 4,194 genes (expression of 2,313 was upregulated and expression of 1,881 was downregulated; Supplementary Table 1), 446 IRGs (expression of 387 was upregulated and expression of 59 was downregulated; Supplementary Table 2) and 36 TFs (expression of 29 was upregulated and expression of seven was downregulated; Table 2) were identified as differentially expressed in PAAD tissues (set threshold: P < 0.01, fold change > 2). The results mentioned above were shown as a heatmap and volcano map (Figure 2).

TABLE 2

Table 2. Differentially expressed transcription factors (TFs).

FIGURE 2

Figure 2. Identification of differentially expressed genes (DEGs), immune-related genes (IRGs), and transcription factors (TFs) in Pancreatic adenocarcinoma (PAAD) vs. normal tissues. (A) Volcano plot revealing clusters of DEGs with upregulated and downregulated expression. (B) The distinction between DEG expression in tumor tissues and normal tissues revealed by a hierarchical clustering heatmap. (C) Volcano plot demonstrates clusters of differentially expressed immune-related genes (DEIRGs) with upregulated and downregulated expression. (D) Heatmap showing the distinction between expression of DEIRGs in tumor tissues and normal tissues. (E) Volcano plot showing clusters of differentially expressed transcription factors (DETFs) with upregulated and downregulated expression. (F) Discrimination between DETFs expression in tumor tissues and normal tissues revealed by a heatmap.

Functional and Pathway Analyses Using GO and KEGG Databases

We wished to elucidate the biological properties and pathways of DEIRGs in PAAD patients. Hence, the GO and KEGG databases were employed. Inevitably, the DEIRGs were enriched in several immune-related molecular functions. The correlation between the top-five most important GO terms and their related DEIRGs was displayed (adjusted-P < 0.05; Figures 3A–C). Among them, “GO: 0019814 immunoglobulin complex” was the most prominent GO term. Figure 3D displays the top-20 significant pathways. The “pathway-DEIRGs” network (Figure 3E) was used for visualizing the reciprocity between the top-10 significant pathways and DEIRGs. Supplementary Table 4 shows 57 significant pathways according to the KEGG database. Adjusted-P < 0.05 was considered indicative of significance. Based on visualized data mining, hsa04060 (“cytokine–cytokine receptor interaction”), hsa04061 (“viral protein interaction with cytokine and cytokine receptor”), and hsa04062 (“chemokine signaling pathway”) were used more often to validate our findings using the KEGG database.

FIGURE 3

Figure 3. Functional-enrichment analyses of DEIRGs in PAAD. (A) The outer circle shows expression (log FC) of DEIRGs in each enriched Gene Ontology (GO) term: red dots which are on each GO term indicate upregulation of DEIRGs. Blue dots indicate downregulated DEIRGs. The inner circle shows the prominence of GO terms (log10-adjusted P-values). (B) The circle represents the relationship between the top-five most significant GO terms and their related DEIRGs. (C) The top-five most significant GO terms and their annotations. (D) The top-20 pathways enriched in DEIRGs are shown in the bubble plot. (E) The top-10 pathways and the corresponding DEIRGs. The blue rectangles represent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The red ellipses indicate upregulated DEIRGs. The green ellipses indicate downregulated DEIRGs.

Regulatory Network of TFs

Univariate Cox regression analysis revealed that 43 DEIRGs were associated with overall survival (OS) (P < 0.01): 37 high-risk IRGs and six low-risk IRGs (Figure 4A). We constructed a regulatory network based on 54 DEIRGs and 36 DETFs (set threshold: P < 0.001; correlation coefficient > 0.4). According to the cutoff criteria, 29 prognostic-related DEIRGs and 14 DETFs (Figure 4B) participated in the establishment of the network. Finally, the regulatory network was constructed and visualized using Cytoscape software¹⁴. The TFs lamin B1 (LMNB1) and lymphoid enhancer-binding factor 1 (LEF1) act as negative regulators of IRG SHC adaptor protein 2 (SHC2) (Figure 4B). Besides, the TF vitamin-D receptor (VDR) had a negative relationship with IRG fibroblast growth factor 17 (FGF17) and neuregulin 2 (NRG2). The specific regulatory relationship between TFs and OS-related IRGs in PAAD was listed in Supplementary Table 5.

FIGURE 4

Figure 4. Overall survival (OS)-related DEIRGs and TFs-IRGs regulatory network. (A) The forest plot of OS-related DEIRGs in PAAD. Red and green dots indicate high risk and low risk, respectively. (B) Regulatory network between prognosis-related DEIRGs and DETFs in PAAD. The red and green circles indicate high-risk and low-risk DEIRGs, respectively. The blue triangles indicate DETFs. The red and green lines represent positive and negative correlation, respectively.

Construction of an Eight-IRGs Prognostic Model

Least absolute shrinkage and selection operator Cox regression analysis was applied to build a prognostic model based on the expression profile of the 43 prognostic DEIRGs mentioned above. Finally, eight genes were selected to construct a prognostic model based on the optimal value of λ (Supplementary Figure 1). The specific formula for the calculation was:

$\text{Risk}\text{score}={\text{e}}^{\begin{array}{c}[(-0.1301)\times \mathrm{expression}\mathrm{of}\mathrm{WFIKKN1}\hfill \\ +0.0016\times \mathrm{expression}\mathrm{of}\mathrm{PLAU}+0.0004\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{OASL}+\left(-0.2278\right)\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{FGF17}+\left(-0.3925\right)\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{NPPA}+0.0258\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{IL20}\text{RB}+0.0251\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{MET}+(-0.0081)\hfill \\ \times \mathrm{expression}\mathrm{of}\mathrm{SHC2}]\hfill \end{array}}$

Pancreatic adenocarcinoma patients were separated into a high-risk group (n = 81) and a low-risk group (n = 81) based on the median value of the risk score (Figure 5C). PCA was undertaken to study the differences between low- and high-risk populations using the expression profiles of all genes, IRGs, and risk-related genes (Figure 6). We discovered that low- and high-risk groups were distributed in different directions (Figure 6C). Patients with high risk were more likely to die sooner than those with a low risk (Figure 5D). The Kaplan–Meier curve demonstrated that patients with high risk showed markedly worse OS than those with low risk (P < 0.001; Figure 5A). The area under the time-dependent ROC curves for 1-, 2-, and 3-years OS reached 0.750, 0.697, and 0.707 respectively. Hence, the predictive performance of the prognostic model exhibited good sensitivity and specificity (Figure 5B). Also, the immune-based prognostic model was relatively consistent. Figure 5E shows the expression of eight IRGs in the form of a heatmap.

FIGURE 5

Figure 5. Prognostic value of eight DEIRGs in PAAD patients. (A) Analyses of Kaplan–Meier curves for OS in PAAD patients using the signature of eight DEIRGs. (B) Receiver operating characteristic (ROC) curve suggesting the feasibility of our prognostic model. (C) Patients in high-risk (red dots) and low-risk (green dots) groups and the distribution of their corresponding risk score. (D) Patients in high-risk (red dots) and low-risk (green dots) groups, and their corresponding survival status. (E) Discrimination of expression of eight prognosis-related IRGs between high-risk and low-risk groups as revealed by a heatmap.

FIGURE 6

Figure 6. Principal component analysis between low-risk and high-risk groups based on different classification methods. (A) All genes, (B) Immune genes, and (C) Risk genes.

Independent Prognostic Value of the Eight-Gene Signature

The independent predictive value of the prognostic signature was assessed by univariate and multivariate Cox regression analyses. Univariate Cox prognostic analyses demonstrated that the risk score was correlated significantly with OS [hazard ratio (HR) = 4.910, 95% confidence interval (CI) = 3.021–7.980, P < 0.001] (Figure 7A). After the multivariate analysis, the risk score remained an independent prognostic factor correlated with OS (HR = 4.868, 95% CI = 2.899–8.175, P < 0.001; Figure 7B). Moreover, univariate and multivariate independent prognostic analyses (Figures 7A,B) showed that the residual tumor and outcomes were also significant independent prognostic factors for survival (P < 0.05; Table 3).

FIGURE 7

Figure 7. Univariate and multivariate independent prognostic analysis in PAAD. (A,B) Forest plots of univariate and multivariate independent prognostic analysis.

TABLE 3

Table 3. Univariate and multivariate independent prognostic analyses.

Eight-IRG Prognostic Model and Clinical Characteristics

Relationships between eight IRGs in the risk-score model and clinical features (age, gender, pathological TNM stage, histology grade, residual tumor, and outcomes) were assessed via the beeswarm packages in R (P < 0.05; Table 4). The cutoff value was determined by the median of the expression of the selected genes. As observed from Figure 8, the median values in the age ≤ 65 group were higher than those in the age > 65 group between mesenchymal epithelial transition factor (MET) expression and riskscore (Figures 8A,B). The median value of the SHC2 and interleukin 20 receptor subunit beta (IL20RB) expression in pathological stage I-II was higher than that in stage III-IV (Figures 8C,D). With regard to histology grade, the median value of MET expression and riskscore in grades 1 and 2 was lower than that in grade 3 and 4, and the trend was exactly the opposite for Fibroblast Growth Factor (FGF)17 (Figures 8E–G). The median values in T1–2 staging were lower than those in T3–4 staging among 2′–5′-oligoadenylate synthetase like (OASL) expression, MET expression, and riskscore (Figures 8H–J). Moreover, the median value of MET expression was lower in residual tumor R0 than that in R1 and 2 (Figure 8K). Additionally, the median values of SHC2, plasminogen activator, urokinase (PLAU) expression, MET expression, IL20RB expression, and riskscore were notably different in PAAD at outcomes CR relative to those at PR + PD + SD (P < 0.05, Figures 8L–P).

TABLE 4

Table 4. Correlation between clinical features.

FIGURE 8

Figure 8. Relationships between the clinical-pathological characteristics and expression of DEIRGs in PAAD. (A) Differences in expression of DEIRGs between the pathological TNM stages I and II/III and IV in PAAD. (B,C) Differences in expression of DEIRGs between the histology G1 and 2/G3 and 4 grades in PAAD. (D,E) Differences in expression of DEIRGs between the T stages T1 and 2/T3 and 4 in PAAD. (F–K) Differences in expression of DEIRGs between the M stages M0/M1/MX in PAAD. (L) Differences in expression of DEIRGs between the residual tumor R0/R1 and R2 in PAAD. (M–P) Differences in expression of DEIRGs between the outcomes (CR/PR + PD + SD) in PAAD.

Immunocyte Infiltration

We wished to ascertain if the eight-IRG prognostic model reflected the status of the PAAD immune microenvironment precisely. Hence, correlation analysis was done to explore the relationship between prognostic IRGs and infiltration of immune cells (Figure 9). The number of DCs, neutrophils, and CD8⁺ T cells was positively correlated with the risk-score prediction model (P < 0.05; Figures 9C,D,F) but the trend of CD4⁺ T cells was opposite (Figure 9B).

FIGURE 9

Figure 9. Relationships between prognostic value and degree of infiltration of six types of immune cells. The relationship of the eight-IRG prognostic model with (A) B cells, (B) CD4 T cells, (C) CD8 T cells, (D) dendritic cells, (E) macrophages, and (F) neutrophils is revealed by scatter diagrams.

External Verification of the Eight-IRG Prognostic Model

Out of the eight prognostic IRGs in our model, the expression of four IRGs was upregulated and that of the remaining four IRGs was downregulated in the TCGA–PAAD cohort. In addition, a GEO dataset (GSE71729) was used to externally verify the difference in expression of eight IRGs between tumor tissues and normal tissues. As expected, the expression of IL20RB, MET, OASL, and PLAU in tumor tissues was significantly higher than that in normal tissues. FGF17, natriuretic peptide A (NPPA), SHC2, and WAP, follistatin/kazal, immunoglobulin, kunitz, and netrin domain containing 1 (WFIKKN1) was not expressed or at a minimal level in tumor tissues (Figure 10A). The protein distribution and expression of FGF17, MET, and SHC2 are displayed in Figures 11A–F, whereas the other five IRGs remained inaccessible in the Human Protein Atlas. Verification using the GEO database further confirmed that PAAD patients in the low-risk group showed a significant OS benefit compared with that in PAAD patients in the high-risk group. The Kaplan–Meier estimator effectively distinguished different groups of various risk (P < 0.01; Figure 10B). The predictive capacity of the signature was confirmed by analyses of the ROC curve. Our results showed that the prognostic signatures of the GEO dataset also performed well in forecasting 1-, 2-, and 3-years survival (Figure 10C).

FIGURE 10

Figure 10. Validation of the eight-gene signature in the Gene Expression Omnibus (GEO) database. (A) Expression of eight IRGs in tumor tissues and normal tissues in the GSE71729 database. (B) Kaplan–Meier curves for low- and high-risk groups in the GSE71729 database (P < 0.01). (C) ROC curve for predicting survival from PAAD based on the risk score of the GSE71729 database. *p < 0.05; **p < 0.01; ***p < 0.001.

FIGURE 11

Figure 11. Representative immunohistochemistry images for expression of FGF17, MET, and SHC2 in pancreatic cancer tissues and normal tissues were shown with the fraction of samples with antibody staining/protein expression level high, medium, low, or not detected. (A–B) Expression of FGF17 in PAAD tissues were lower than that in normal tissues. (C–D) Expression of MET in PAAD tissues were higher than that in normal tissues. (E–F) Expression of SHC2 in PAAD tissues were obviously lower than that in normal tissues.

Discussion

Pancreatic cancer remains a lethal type of cancer due to its poor prognosis and lack of efficacious therapeutic approaches. Precise prediction of OS after contracting PAAD is very important for the choice of therapeutic method and improving the prognosis.

Pancreatic cancer lacks reliable and effective prognostic biomarkers related to the tumor immune microenvironment. An effective prediction model to accurately assess the prognosis of PAAD is long overdue. We intended to explore DEIRGs and establish a model of PAAD based on IRGs to uncover the biomarkers that predict the diagnosis and prognosis of PAAD.

In our study, 446 DEIRGs of PC were identified by comprehensive analyses. Analyses of pathway enrichment revealed that these DEIRGs correlated with the inflammatory response and typical tumor-related pathways shown in Supplementary Table 4. Most of them were related to the progression and treatment of PAAD. Cytokines and their correlated pathways may play a relevant part in PAAD progression and immune evasion (Padoan et al., 2019; Dey et al., 2020). As a vital component of the signaling between cancer cells and surrounding stromal cells, chemokine signaling participates in the development of the supportive TME of PAAD (Sleightholm et al., 2017). The Janus kinase family/signal transducer and activator of transcription (JAK/STAT) signaling pathway were central to tumor growth, tumor survival, and systemic inflammation, particularly in PC (Quintás-Cardama and Verstovsek, 2013; von Ahrens et al., 2017). In addition, two studies (Hurwitz et al., 2015; Beatty et al., 2019) showed that inhibitors of the JAK/STAT pathway may have clinical benefit. A follow-up study of a general population indicated that the high cytotoxic activity of natural killer (NK) cells is linked to a reduced risk of cancer (Imai et al., 2000). Lee et al. (2020) stated that the activity of NK cells decreased as cancer progressed, and that decreased activity of NK cells was associated with poor clinical outcomes. NF-κB is a pro-inflammatory signaling pathway in pancreatitis and PAAD. Increased basal levels and/or inducible levels of NF-κB activation are strongly linked to several aspects of treatment resistance, as well as the proliferation and metastasis of tumor cells in PAAD (Arlt et al., 2012; Kabacaoglu et al., 2019). In addition, NF-κB-mediated chemokine signaling plays a crucial part in the therapy resistance of PC (Geismann et al., 2019). Signaling by T-helper (Th1) and Th2 cytokines is complex in the microenvironment of pancreatic tumors (Andrianifahanana et al., 2006). The presence of tumor-infiltrating lymphocytes with high Th2:Th1 ratios demonstrates a poor prognosis in PAAD (Seicean et al., 2009). He et al. (2011) revealed that the accumulation of Th17 cells and their relevant cytokine levels in PC tissues may manifest engagement in the invasion and metastasis of PC, which may thereby have an impact on the prognosis. We conducted a comprehensive investigation of the biological functions of DEIRGs in PAAD populations to provide a basis for elucidating their possible molecular regulatory mechanisms.

More and more studies have found that abnormally expressed TFs in tumor tissues were related to aggressive diseases and poor prognosis. The research on new drugs that target specific TFs had great potential in developing clinically relevant strategies for the treatment of malignant tumors (Sankpal et al., 2012). To further investigate the possible molecular regulatory mechanisms, a TFs-mediated prognosis-related IRGs network was structured to find the significant TFs regulating DEIRGs in this network. DETFs such as basic helix-loop-helix family member E40 (BHLHE40), E2F transcription factor 1 (E2F1), early growth response 2 (EGR2), FOS like 1, AP-1 transcription factor subunit (FOSL1), forkhead box M1 (FOXM1), kruppel like factor 5 (KLF5), LEF1, LMNB1, peroxisome proliferator activated receptor gamma (PPARG), PR/SET domain 1 (PRDM1), transcription factor AP-2 alpha (TFAP2A), tumor protein P63 (TP63), and VDR might regulate the DEIRGs in PAAD. BHLHE40 expression was upregulated by transforming growth factor (TGF)-β, and affected the morphology, migration, and invasion of PC cells by changing the expression of factors related to epithelial-to-mesenchymal (EMT) transition (Wu et al., 2012). E2F1-mediated overexpression of long non-coding (lnc)RNA-pancreatic cancer associated transcript 1 (PLACT1) promotes the growth of PAAD by continuously activating the NF-κB pathway and forming a positive feedback loop with IκBα in PC (Ren et al., 2020). Vallejo et al. (2017) revealed FOSL1 to be an oncogene in KRAS-driven lung cancer and PC, which partially factors through transcriptional regulation of a subset of genes involved in the mitotic machinery. Zhou et al. (2019) revealed an important epigenetic modification to FOXM1, and increased expression of FOXM1 suppressed the maturation of bone marrow−derived DCs via direct activation of Wnt5a signaling pathway and weakened the promotion of T−cell proliferation. He et al. (2018) demonstrated that KLF5 depletion in oncogenic Kras-expressing mouse PC cells reduced proliferation of tumor cells and PC progression. TP63 is a member of the p53 family and is transcribed from two promoters to produce two subtypes: TAp63 and ΔNp63 (Gonfloni et al., 2015). TP63 reprograms enhancers to drive squamous transdifferentiation in PC (Somerville et al., 2020). Sherman and collaborators discovered that the VDR is expressed in the stroma from PC cells and acts as a “master” transcriptional regulator of pancreatic stellate cells, thereby resulting in induced transcriptional reprogramming of tumor stroma in PAAD.

We innovatively established a TFs-mediated prognosis-related IRGs regulatory network in PAAD by bioinformatics analysis. This network showed that TFs regulated IRGs positively and negatively, which supplied a novel method to explore the IRGs underlying regulatory mechanisms in PAAD at the molecular level.

Eight IRGs involved in the prognostic model were considered to be potential biomarkers in PAAD. Among the eight genes, MET, OASL, SHC2, and PLAU have been well studied in PAAD compared with other IRGs. Nan and coworkers found that hepatocyte growth factor (HGF) promotes the invasion and migration of PC cells by activating the HGF/c-Met pathway (Nan et al., 2019). Besides, MET/HGF co-targeting may represent a treatment option for patients with PC (Modica et al., 2018). As a member of the OAS protein family, OASL is associated with the innate immune defense against viral infections. Glaß et al. (2020) identified OASL to be a candidate oncogenic RNA-binding protein with partially validated target potential in PC. Recently, PLAU has been reported to be an oncogene that activates EMT progression in PAAD (Zhao X. et al., 2020). SHC2 was a proverbial adaptor molecule that binds to receptor tyrosine kinases via its SH2 domain. Teodorczyk and coworkers reported that CD95L could induce SCK recruitment and activation of the phosphoinositide 3-kinase/extracellular signal-regulated kinase (PI3K/ERK) pathway by stimulating CD95 receptors and, ultimately, lead to PC cell-cycle progression (Teodorczyk et al., 2015). One review stated that high expression of IL20RB was related to poor survival, thereby suggesting its oncogenic potential in PAAD (Haider et al., 2014). FGF17 was a member of the FGF8 subfamily, which promotes the development and progression of hepatocellular carcinoma (Gauglhofer et al., 2011). In addition, FGF17 was overexpressed in human prostate cancer, and involved in the progression of prostate cancer to high−grade disease (Heer et al., 2004). Both studies have reported that FGF17 may be a novel tumor-promoting gene whose expression is upregulated in neoplasms, data which contradicted our findings. The exact role of FGF17 in PAAD is not known. Few related studies have reported NPPA or WFIKKN1 being involved in PAAD.

Our TFs-IRGs-mediated network contained five of the eight modeling genes: PLAU, OASL, FGF17, MET, and SHC2. Their interactions with tumor-associated TFs can provide a certain theoretical direction/basis for mechanistic studies. Therefore, further study of the potential regulatory mechanisms of these prognostic immune genes in PAAD is needed.

To clarify the immune microenvironment in PAAD, a correlation analysis on immunocytes was done based on the TIMER database. Results indicated that lower infiltration of DCs, CD8⁺ T cells, and neutrophils may be observed in low-risk patients, whereas the tendency of CD4⁺ T cells was the opposite. DCs, neutrophils, and CD8⁺ T cells exhibited a significantly positive regulatory relationship with the prognostic model. Thus, our model may act as a predictive factor for increased infiltration of immune cells. One study reported that higher numbers of CD4⁺ T lymphocytes were significantly associated with longer survival, which echoed our findings (Ino et al., 2013). A recent study showed that intratumoral infiltration by CD8⁺ T lymphocytes and neutrophils and a favorable prognosis in PAAD patients were tightly linked (Miksch et al., 2019), which is the reverse of our results. Thus, our results must be validated by further investigations. Whether the infiltration level of DCs in tumors indicates the clinical prognosis of PAAD patients has not been reported. Studies on other tumors have yielded inconsistent or even conflicting results that doubt the value of infiltrating DCs (Karthaus et al., 2012). The exact role of immunocytes in PAAD has not been clarified. Considering the different levels of immunocyte infiltration between high-risk and low-risk PAAD groups, suitable immunotherapy strategies can be selected based on the basis of the immune microenvironment in PAAD.

So far, several studies have proposed that prognostic gene signatures based on mRNA levels can predict the OS of PC prognosis. For instance, Birnbaum et al. (2017) built a 25-gene classifier that helped select patients with resectable disease for immediate surgery or neoadjuvant chemotherapy. Another study established a four-gene signature for prediction of OS from PC based on gene-expression data from the GEO database [1-, 2-, and 3-years survival area under the curve (AUC) reached 0.715, 0.654, and 0.715, respectively] (Yan et al., 2019). A recent study investigated the survival-associated genes from the integrated analysis of multiple datasets, and established prognostic signatures in PAAD (1-, 2-, and 3-years survival AUC reached 0.699, 0.637, and 0.621, respectively; Wu et al., 2019). At present, there are few studies on the relationship between IRGs and the prognosis of PAAD. The latest research developed an immune prognostic model to identify low-risk patients who may benefit from immunotherapy (Gu et al., 2021). However, this predictive model still lacks an external cohort to verify the effectiveness of the model. We used a specialized immune database to explore the relationship between many IRGs and the prognosis of PAAD patients. Subsequently, we established a new immune-related prognostic signature. No overlap was found between the eight-gene signatures we developed and the one defined previously. Besides, the riskscore had a robust predictive performance with 1-, 2-, and 3-years survival AUC reached 0.750, 0.697, and 0.704, respectively. The predictive performance of our prognostic model was superior or comparable with that reported in other studies, and this model prediction is verified in an external validation cohort. These results suggest that an immune-related prognostic signature may be a valid marker for the prediction of the PC prognosis.

Nevertheless, our study still has perceived limitations. Firstly, we only used data from sections of public databases to build and validate our prediction model. Therefore, one must conduct more prospective studies to verify its clinical applicability. Secondly, we excluded many prominent prognostic genes in PAAD, so the potential weakness inherent in constructing a prognostic model with a single hallmark is inevitable. Moreover, the protein expression of IRGs related to the prognosis, and their potential molecular mechanisms in the pathogenesis and development of PAAD, must be confirmed by additional experimental studies.

Conclusion

We defined a novel eight-IRG model as an independent prognostic predictor for PAAD. The prognostic value of this model was verified by an external validation database. Moreover, the correlation between the eight-IRG prognostic model and infiltrated immunocytes could demonstrate its pivotal role in the PAAD immune microenvironment, which could be utilized as a new prognostic and therapeutic biomarker in PAAD patients.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

MM, HL, and YL wrote the main manuscript text. MM and HL prepared figures and tables. YC, BX, and RZ reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by Joint Funds for the Innovation of Science and Technology, Fujian province (2017Y9061).

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

The authors would like to thank the TCGA, GEO, GTEx, ImmPort, Cistrome Cancer, HPA, and TIMER databases for the availability of the data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.702102/full#supplementary-material

Footnotes

1. ^ www.gtexportal.org/
2. ^ https://www.ncbi.nlm.nih.gov/geo/
3. ^ www.immport.org/
4. ^ www.cistrome.org/
5. ^ www.bioconductor.org/packages/release/bioc/html/limma.html
6. ^ https://david.ncifcrf.gov/
7. ^ https://cran.r-project.org/web/packages/GOplot/citation.html
8. ^ www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html
9. ^ https://cran.r-project.org/web/views/Survival.html
10. ^ https://cran.r-project.org/web/packages/glmnet/index.html
11. ^ https://cran.r-project.org/web/packages/survivalROC/index.html
12. ^ www.proteinatlas.org/
13. ^ https://cistrome.shinyapps.io/timer/
14. ^ www.cytoscape.org/

References

Andrianifahanana, M., Chauhan, S. C., Choudhury, A., Moniaux, N., Brand, R. E., Sasson, A. A., et al. (2006). MUC4-expressing pancreatic adenocarcinomas show elevated levels of both T1 and T2 cytokines: potential pathobiologic implications. Am. J. Gastroenterol. 101, 2319–2329. doi: 10.1111/j.1572-0241.2006.00871.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Arlt, A., Schäfer, H., and Kalthoff, H. (2012). The ‘N-factors’ in pancreatic cancer: functional relevance of NF-κB, NFAT and Nrf2 in pancreatic cancer. Oncogenesis 1:e35. doi: 10.1038/oncsis.2012.35

PubMed Abstract | CrossRef Full Text | Google Scholar

Beatty, G. L., Shahda, S., Beck, T., Uppal, N., Cohen, S. J., Donehower, R., et al. (2019). A Phase Ib/II Study of the JAK1 inhibitor, itacitinib, plus nab-paclitaxel and gemcitabine in advanced solid tumors. Oncologist 24, 14–e10. doi: 10.1634/theoncologist.2017-0665

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Andorf, S., Gomes, L., Dunn, P., Schaefer, H., Pontius, J., et al. (2014). ImmPort: disseminating data to the public for the future of immunology. Immunol. Res. 58, 234–239. doi: 10.1007/s12026-014-8516-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Bianco, A., Perrotta, F., Barra, G., Malapelle, U., Rocco, D., and De Palma, R. (2019). Prognostic factors and biomarkers of responses to immune checkpoint inhibitors in lung cancer. Int. J. Mol. Sci. 20:4931. doi: 10.3390/ijms20194931

PubMed Abstract | CrossRef Full Text | Google Scholar

Birnbaum, D. J., Finetti, P., Lopresti, A., Gilabert, M., Poizat, F., Raoul, J. L., et al. (2017). A 25-gene classifier predicts overall survival in resectable pancreatic cancer. BMC Med. 15:170. doi: 10.1186/s12916-017-0936-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Carithers, L. J., Ardlie, K., Barcus, M., Branton, P. A., Britton, A., Buia, S. A., et al. (2015). A novel approach to high-quality postmortem tissue procurement: the GTEx project. Biopreserv. Biobank. 13, 311–319. doi: 10.1089/bio.2015.0032

PubMed Abstract | CrossRef Full Text | Google Scholar

Dennis, G. Jr., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 4:R60.

Google Scholar

Dey, P., Li, J., Zhang, J., Chaurasiya, S., Strom, A., Wang, H., et al. (2020). Oncogenic KRAS-driven metabolic reprogramming in pancreatic cancer cells utilizes cytokines from the tumor microenvironment. Cancer Discov. 10, 608–625. doi: 10.1158/2159-8290.cd-19-0297

PubMed Abstract | CrossRef Full Text | Google Scholar

Erkan, M., Hausmann, S., Michalski, C. W., Fingerle, A. A., Dobritz, M., Kleeff, J., et al. (2012). The role of stroma in pancreatic cancer: diagnostic and therapeutic implications. Nat. Rev. Gastroenterol. Hepatol. 9, 454–467. doi: 10.1038/nrgastro.2012.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Partensky, C., and Bray, F. (2016). More deaths from pancreatic cancer than breast cancer in the EU by 2017. Acta Oncol. (Stockholm, Sweden) 55, 1158–1160. doi: 10.1080/0284186x.2016.1197419

PubMed Abstract | CrossRef Full Text | Google Scholar

Gauglhofer, C., Sagmeister, S., Schrottmaier, W., Fischer, C., Rodgarkia-Dara, C., Mohr, T., et al. (2011). Up-regulation of the fibroblast growth factor 8 subfamily in human hepatocellular carcinoma for cell survival and neoangiogenesis. Hepatology (Baltimore, Md) 53, 854–864. doi: 10.1002/hep.24099

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, P., Wang, W., Li, L., Zhang, G., Gao, Z., Tang, Z., et al. (2019). Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of colorectal cancer. Biomed. Pharmacother. = Biomed. Pharmacother. 118:109228. doi: 10.1016/j.biopha.2019.109228

PubMed Abstract | CrossRef Full Text | Google Scholar

Geismann, C., Schäfer, H., Gundlach, J. P., Hauser, C., Egberts, J. H., Schneider, G., et al. (2019). NF-κB dependent chemokine signaling in pancreatic cancer. Cancers 11:1445. doi: 10.3390/cancers11101445

PubMed Abstract | CrossRef Full Text | Google Scholar

Glaß, M., Michl, P., and Hüttelmaier, A. S. (2020). RNA Binding proteins as drivers and therapeutic target candidates in pancreatic ductal adenocarcinoma. Int. J. Mol. Sci. 21:4190. doi: 10.3390/ijms21114190

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonfloni, S., Caputo, V., and Iannizzotto, V. (2015). P63 in health and cancer. Int. J. Dev. Biol. 59, 87–93. doi: 10.1387/ijdb.150045sg

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, X., Zhang, Q., Wu, X., Fan, Y., and Qian, J. (2021). Gene coexpression network approach to develop an immune prognostic model for pancreatic adenocarcinoma. World J. Surgical Oncol. 19:112. doi: 10.1186/s12957-021-02201-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Haider, S., Wang, J., Nagano, A., Desai, A., Arumugam, P., Dumartin, L., et al. (2014). A multi-gene signature predicts outcome in patients with pancreatic ductal adenocarcinoma. Genome Med. 6:105. doi: 10.1186/s13073-014-0105-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanley, J. A., and McNeil, B. J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 143, 29–36. doi: 10.1148/radiology.143.1.7063747

PubMed Abstract | CrossRef Full Text | Google Scholar

He, P., Yang, J. W., Yang, V. W., and Bialkowska, A. B. (2018). Krüppel-like Factor 5, increased in pancreatic ductal adenocarcinoma, promotes proliferation, acinar-to-ductal metaplasia, pancreatic intraepithelial neoplasia, and tumor growth in mice. Gastroenterology 154, 1494.e13–1508.e13. doi: 10.1053/j.gastro.2017.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

He, S., Fei, M., Wu, Y., Zheng, D., Wan, D., Wang, L., et al. (2011). Distribution and clinical significance of Th17 cells in the tumor microenvironment and peripheral blood of pancreatic cancer patients. Int. J. Mol. Sci. 12, 7424–7437. doi: 10.3390/ijms12117424

PubMed Abstract | CrossRef Full Text | Google Scholar

Heer, R., Douglas, D., Mathers, M. E., Robson, C. N., and Leung, H. Y. (2004). Fibroblast growth factor 17 is over-expressed in human prostate cancer. J. Pathol. 204, 578–586. doi: 10.1002/path.1668

PubMed Abstract | CrossRef Full Text | Google Scholar

Higuera, O., Ghanem, I., Nasimi, R., Prieto, I., Koren, L., and Feliu, J. (2016). Management of pancreatic cancer in the elderly. World J. Gastroenterol. 22, 764–775. doi: 10.3748/wjg.v22.i2.764

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, R., Mao, M., Lu, Y., Yu, Q., and Liao, L. (2020). A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging 12, 6966–6980. doi: 10.18632/aging.103054

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurwitz, H. I., Uppal, N., Wagner, S. A., Bendell, J. C., Beck, J. T., and Wade, S. M. 3rd, et al. (2015). Randomized, double-blind, phase II study of ruxolitinib or placebo in combination with capecitabine in patients with metastatic pancreatic cancer for whom therapy with gemcitabine has failed. J. Clin. Oncol. : Off. J. Am. Soc. Clin. Oncol. 33, 4039–4047. doi: 10.1200/jco.2015.61.4578

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilic, M., and Ilic, I. (2016). Epidemiology of pancreatic cancer. World J. Gastroenterol. 22, 9694–9705. doi: 10.3748/wjg.v22.i44.9694

PubMed Abstract | CrossRef Full Text | Google Scholar

Imai, K., Matsuyama, S., Miyake, S., Suga, K., and Nakachi, K. (2000). Natural cytotoxic activity of peripheral-blood lymphocytes and cancer incidence: an 11-year follow-up study of a general population. Lancet (London, England) 356, 1795–1799. doi: 10.1016/s0140-6736(00)03231-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Incio, J., Liu, H., Suboj, P., Chin, S. M., Chen, I. X., Pinter, M., et al. (2016). Obesity-induced inflammation and desmoplasia promote pancreatic cancer progression and resistance to chemotherapy. Cancer Discov. 6, 852–869. doi: 10.1158/2159-8290.cd-15-1177

PubMed Abstract | CrossRef Full Text | Google Scholar

Ino, Y., Yamazaki-Itoh, R., Shimada, K., Iwasaki, M., Kosuge, T., Kanai, Y., et al. (2013). Immune cell infiltration as an indicator of the immune microenvironment of pancreatic cancer. Br. J. Cancer 108, 914–923. doi: 10.1038/bjc.2013.32

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Xie, J., Huang, W., Chen, H., Xi, S., Han, Z., et al. (2019). Tumor immune microenvironment and chemosensitivity signature for predicting response to chemotherapy in gastric cancer. Cancer Immunol. Res. 7, 2065–2073. doi: 10.1158/2326-6066.cir-19-0311

PubMed Abstract | CrossRef Full Text | Google Scholar

Kabacaoglu, D., Ruess, D. A., Ai, J., and Algül, H. (2019). NF-κB/Rel transcription factors in pancreatic cancer: focusing on RelA, c-Rel, and RelB. Cancers 11:937. doi: 10.3390/cancers11070937

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361. doi: 10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

Karthaus, N., Torensma, R., and Tel, J. (2012). Deciphering the message broadcast by tumor-infiltrating dendritic cells. Am. J. Pathol. 181, 733–742. doi: 10.1016/j.ajpath.2012.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, Y., Feng, Z. C., Zhang, Y. L., Liu, X. F., Ma, Y., Zhao, Z. M., et al. (2020). Identification of immune-related genes contributing to the development of glioblastoma using weighted gene co-expression network analysis. Front. Immunol. 11:1281. doi: 10.3389/fimmu.2020.01281

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H. S., Leem, G., Kang, H., Jo, J. H., Chung, M. J., Jang, S. J., et al. (2020). Peripheral natural killer cell activity is associated with poor clinical outcomes in pancreatic ductal adenocarcinoma. J. Gastroenterol. Hepatol. 36, 516–522. doi: 10.1111/jgh.15265

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77, e108–e110. doi: 10.1158/0008-5472.can-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, C., Shi, S., Meng, Q., Liang, D., Ji, S., Zhang, B., et al. (2017). Complex roles of the stroma in the intrinsic resistance to gemcitabine in pancreatic cancer: where we are and where we are going. Exp. Mol. Med. 49:e406. doi: 10.1038/emm.2017.255

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S. L., Bian, L. J., Liu, Z. X., Chen, Q. Y., Sun, X. S., Sun, R., et al. (2020). Development and validation of the immune signature to predict distant metastasis in patients with nasopharyngeal carcinoma. J. Immunother. Cancer 8:e000205. doi: 10.1136/jitc-2019-000205

PubMed Abstract | CrossRef Full Text | Google Scholar

Markowitz, J., Brooks, T. R., Duggan, M. C., Paul, B. K., Pan, X., Wei, L., et al. (2015). Patients with pancreatic adenocarcinoma exhibit elevated levels of myeloid-derived suppressor cells upon progression of disease. Cancer Immunol. Immunother.: CII 64, 149–159. doi: 10.1007/s00262-014-1618-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, S., Meyer, C. A., Zheng, R., Qin, Q., Wu, Q., Jiang, P., et al. (2017). Cistrome cancer: a web resource for integrative gene regulation modeling in cancer. Cancer Res. 77, e19–e22. doi: 10.1158/0008-5472.can-17-0327

PubMed Abstract | CrossRef Full Text | Google Scholar

Miksch, R. C., Schoenberg, M. B., Weniger, M., Bösch, F., Ormanns, S., Mayer, B., et al. (2019). Prognostic impact of tumor-infiltrating lymphocytes and neutrophils on survival of patients with upfront resection of pancreatic cancer. Cancers 11:39. doi: 10.3390/cancers11010039

PubMed Abstract | CrossRef Full Text | Google Scholar

Modica, C., Tortarolo, D., Comoglio, P. M., Basilico, C., and Vigna, E. M. E. T. (2018). /HGF Co-targeting in pancreatic cancer: a tool to provide insight into the tumor/stroma crosstalk. Int. J. Mol. Sci. 19:3920. doi: 10.3390/ijms19123920

PubMed Abstract | CrossRef Full Text | Google Scholar

Nan, L., Qin, T., Xiao, Y., Qian, W., Li, J., Wang, Z., et al. (2019). Pancreatic stellate cells facilitate perineural invasion of pancreatic cancer via HGF/c-Met pathway. Cell Transplantation 28, 1289–1298. doi: 10.1177/0963689719851772

PubMed Abstract | CrossRef Full Text | Google Scholar

Padoan, A., Plebani, M., and Basso, D. (2019). Inflammation and pancreatic cancer: focus on metabolism, cytokines, and immunity. Int. J. Mol. Sci. 20:676. doi: 10.3390/ijms20030676

PubMed Abstract | CrossRef Full Text | Google Scholar

Pontén, F., Schwenk, J. M., Asplund, A., and Edqvist, P. H. (2011). The Human Protein Atlas as a proteomic resource for biomarker discovery. J. Internal Med. 270, 428–446. doi: 10.1111/j.1365-2796.2011.02427.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Quintás-Cardama, A., and Verstovsek, S. (2013). Molecular pathways: jak/STAT pathway: mutations, inhibitors, and resistance. Clin. Cancer Res. : Off. J. Am. Assoc. Cancer Res. 19, 1933–1940. doi: 10.1158/1078-0432.ccr-12-0284

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahib, L., Smith, B. D., Aizenberg, R., Rosenzweig, A. B., Fleshman, J. M., and Matrisian, L. M. (2014). Projecting cancer incidence and deaths to 2030: the unexpected burden of thyroid, liver, and pancreas cancers in the United States. Cancer Res. 74, 2913–2921. doi: 10.1158/0008-5472.can-14-0155

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, B., Cui, M., Yang, G., Wang, H., Feng, M., You, L., et al. (2018). Tumor microenvironment participates in metastasis of pancreatic cancer. Mol. Cancer 17:108. doi: 10.1186/s12943-018-0858-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, X., Chen, C., Luo, Y., Liu, M., Li, Y., Zheng, S., et al. (2020). lncRNA-PLACT1 sustains activation of NF-κB pathway through a positive feedback loop with IκBα/E2F1 axis in pancreatic cancer. Mol. Cancer 19:35. doi: 10.1186/s12943-020-01153-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sankpal, U. T., Maliakal, P., Bose, D., Kayaleh, O., Buchholz, D., and Basha, R. (2012). Expression of specificity protein transcription factors in pancreatic cancer and their association in prognosis and therapy. Curr. Med. Chem. 19, 3779–3786. doi: 10.2174/092986712801661077

PubMed Abstract | CrossRef Full Text | Google Scholar

Seicean, A., Popa, D., Mocan, T., Cristea, V., and Berindan-Neagoe, I. (2009). Th1 and Th2 profiles in patients with pancreatic cancer compared with chronic pancreatitis. Pancreas 38, 594–595. doi: 10.1097/MPA.0b013e31819313d0

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization paths for Cox’s proportional hazards model via coordinate descent. J. Stat. Softw. 39, 1–13. doi: 10.18637/jss.v039.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Sleightholm, R. L., Neilsen, B. K., Li, J., Steele, M. M., Singh, R. K., Hollingsworth, M. A., et al. (2017). Emerging roles of the CXCL12/CXCR4 axis in pancreatic cancer progression and therapy. Pharmacol. Therapeutics 179, 158–170. doi: 10.1016/j.pharmthera.2017.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Somerville, T. D., Biffi, G., Daßler-Plenker, J., Hur, S. K., He, X. Y., Vance, K. E., et al. (2020). Squamous trans-differentiation of pancreatic cancer cells promotes stromal inflammation. eLife 9:e53381. doi: 10.7554/eLife.53381

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Y., Xu, X., Guo, S., Zhang, C., Tang, Y., Tian, Y., et al. (2014). An increased abundance of tumor-infiltrating regulatory T cells is correlated with the progression and prognosis of pancreatic ductal adenocarcinoma. PloS one 9:e91551. doi: 10.1371/journal.pone.0091551

PubMed Abstract | CrossRef Full Text | Google Scholar

Taube, J. M., Galon, J., Sholl, L. M., Rodig, S. J., Cottrell, T. R., Giraldo, N. A., et al. (2018). Implications of the tumor immune microenvironment for staging and therapeutics. Modern Pathol. : Off. J. U.S. Can. Acad. Pathol., Inc 31, 214–234. doi: 10.1038/modpathol.2017.156

PubMed Abstract | CrossRef Full Text | Google Scholar

Teodorczyk, M., Kleber, S., Wollny, D., Sefrin, J. P., Aykut, B., Mateos, A., et al. (2015). CD95 promotes metastatic spread via Sck in pancreatic ductal adenocarcinoma. Cell Death Differ. 22, 1192–1202. doi: 10.1038/cdd.2014.217

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat. Med. 16, 385–395. doi: 10.1002/(sici)1097-0258(19970228)16:4

CrossRef Full Text | Google Scholar

Vallejo, A., Perurena, N., Guruceaga, E., Mazur, P. K., Martinez-Canarias, S., Zandueta, C., et al. (2017). An integrative approach unveils FOSL1 as an oncogene vulnerability in KRAS-driven lung and pancreatic cancer. Nat. Commun. 8:14294. doi: 10.1038/ncomms14294

PubMed Abstract | CrossRef Full Text | Google Scholar

von Ahrens, D., Bhagat, T. D., Nagrath, D., Maitra, A., and Verma, A. (2017). The role of stromal cancer-associated fibroblasts in pancreatic cancer. J. Hematol. Oncol. 10:76. doi: 10.1186/s13045-017-0448-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Walter, W., Sánchez-Cabo, F., and Ricote, M. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics (Oxford, England) 31, 2912–2914. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Lang, M., Zhao, T., Feng, X., Zheng, C., Huang, C., et al. (2017). Cancer-FOXP3 directly activated CCL5 to recruit FOXP3(+)Treg cells in pancreatic ductal adenocarcinoma. Oncogene 36, 3048–3058. doi: 10.1038/onc.2016.458

PubMed Abstract | CrossRef Full Text | Google Scholar

Whatcott, C. J., Diep, C. H., Jiang, P., Watanabe, A., LoBello, J., Sima, C., et al. (2015). Desmoplasia in primary tumors and metastatic lesions of pancreatic cancer. Clin. Cancer Res. : Off. J. Am. Assoc. Cancer Res. 21, 3561–3568. doi: 10.1158/1078-0432.ccr-14-1051

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, M., Li, X., Zhang, T., Liu, Z., and Zhao, Y. (2019). Identification of a nine-gene signature and establishment of a prognostic nomogram predicting overall survival of pancreatic cancer. Front. Oncol. 9:996. doi: 10.3389/fonc.2019.00996

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Sato, F., Yamada, T., Bhawal, U. K., Kawamoto, T., Fujimoto, K., et al. (2012). The BHLH transcription factor DEC1 plays an important role in the epithelial-mesenchymal transition of pancreatic cancer. Int. J. Oncol. 41, 1337–1346. doi: 10.3892/ijo.2012.1559

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, X., Wan, H., Hao, X., Lan, T., Li, W., Xu, L., et al. (2019). Importance of gene expression signatures in pancreatic cancer prognosis and the establishment of a prediction model. Cancer Manag. Res. 11, 273–283. doi: 10.2147/cmar.s185205

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Liu, T., Nan, H., Wang, Y., Chen, H., Zhang, X., et al. (2020). Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma. J. Cell. Physiol. 235, 1025–1035. doi: 10.1002/jcp.29018

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B., Wang, Y., Wang, Y., Chen, W., Liu, P. H., Kong, Z., et al. (2020). Systematic identification, development, and validation of prognostic biomarkers involving the tumor-immune microenvironment for glioblastoma. J. Cell. Physiol. 236, 507–522. doi: 10.1002/jcp.29878

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Liu, Z., Ren, Z., Wang, H., Wang, Z., Zhai, J., et al. (2020). Triptolide inhibits pancreatic cancer cell proliferation and migration via down-regulating PLAU based on network pharmacology of Tripterygium wilfordii Hook F. Eur. J. Pharmacol. 880:173225. doi: 10.1016/j.ejphar.2020.173225

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Z., Chen, H., Xie, R., Wang, H., Li, S., Xu, Q., et al. (2019). Epigenetically modulated FOXM1 suppresses dendritic cell maturation in pancreatic cancer and colon cancer. Mol. Oncol. 13, 873–893. doi: 10.1002/1878-0261.12443

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwing, N., Failmezger, H., Ooi, C. H., Hibar, D. P., Cañamero, M., Gomes, B., et al. (2020). Analysis of spatial organization of suppressive myeloid cells and effector T cells in colorectal cancer-a potential tool for discovering prognostic biomarkers in clinical research. Front. Immunol. 11:550250. doi: 10.3389/fimmu.2020.550250

PubMed Abstract | CrossRef Full Text | Google Scholar