Three Immune-Related Prognostic mRNAs as Therapeutic Targets for Pancreatic Cancer

Objective: Pancreatic cancer is a highly lethal malignancy globally. This study aimed to probe and validate immune-related prognostic mRNAs as therapeutic targets for pancreatic cancer. Methods: Gene transcriptome data of pancreatic cancer and normal pancreas were retrieved from TCGA-GTEx projects. Two thousand four hundred and ninety-eight immune-related genes were obtained from the IMMUPORT database. Abnormally expressed immune-related genes were then identified. Under univariate and multivariate cox models, a gene signature was constructed. Its predictive efficacy was assessed via ROCs. The interactions between the 21 genes were analyzed by Spearson analysis and PPI network. Using the GEPIA and The Human Protein Atlas databases, their expression and prognostic value were evaluated. The TIMER database was utilized to determine the relationships between MET, OAS1, and OASL mRNAs and immune infiltrates. Finally, their mRNA expression was externally verified in the GSE15471 and GSE62452 datasets. Results: An immune-related 21-gene signature was developed for predicting patients' prognosis. Following verification, this signature exhibited the well predictive performance. There were physical and functional interactions between them. MET, OAS1, and OASL mRNAs were all up-regulated in pancreatic cancer and associated with unfavorable prognosis. They showed strong correlations with tumor progression. Furthermore, the three mRNAs were distinctly associated with immune infiltrates. Their up-regulation was confirmed in the two external datasets. Conclusion: These findings identified three immune-related prognostic mRNAs MET, OAS1, and OASL, which may assist clinicians to choose targets for immunotherapy and make personalized treatment strategy for pancreatic cancer patients.


INTRODUCTION
Pancreatic cancer is a highly lethal disease with fairly high mortality globally (1). Annual mortality is nearly equal to incidence in many countries such as China (2). The 5-year survival is ∼7% (3). Its unfavorable prognosis is mainly attributed to local infiltration and distant metastases. Pancreatic ductal adenocarcinoma accounts for around 95% of all cases. Only 10% of cases are due to genetic factors (4). Smoking, drinking, and obesity are common modifiable risk factors for this disease (5). Routine therapy methods like surgery offer dissatisfactory clinical outcomes. Only ∼20% could benefit from curative surgical resection (6). Recently, immunotherapy has become a promising adjunct treatment for pancreatic cancer (7). Nevertheless, most of pancreatic cancer patients are resistant to most therapies including immunotherapy because tumors may evade immune surveillance (8). Moreover, there is currently no targeted therapy against driver genes in pancreatic cancer. Thus, it is of necessity to develop novel strategies to build up treatment efficacy.
High-throughput sequencing may assist us to probe therapeutic targets for cancer therapies and understand the underlying mechanisms of the anti-cancer efficacy in depth (9). Pancreatic cancer is featured by distinct immune disorders. Components in immune system contribute to the initiation and development of pancreatic cancer (10). Integrated analysis of relationships between immune-related genes and clinical outcomes of pancreatic cancer is critical to explore novel prognostic markers as well as therapeutic targets. For example, Wu et al. found that three immune-related genes CKLF, ERAP2, and EREG showed distinct correlations with pancreatic cancer patients' survival (11). In this study, we identified three immune-related prognostic genes MET, OAS1, and OASL from the 21-gene signature. Following multiple dataset verification, these genes could be promising therapeutic targets as well as prognostic markers in pancreatic cancer.

Downloaded Datasets
Gene transcriptome data of 165 normal pancreas samples from the Genotype-Tissue Expression (GTEx) project and 178 pancreatic cancer samples from The Cancer Genome Atlas (TCGA) database were downloaded based on the UCSC Xena browser (https://xenabrowser.net/datapages/). Clinical information of these pancreatic cancer patients including gender, pathologic T, N and histologic grade was retrieved from Genomic Data Commons (GDC). Those without available follow-up data were removed.

Screening for Abnormally Expressed Immune-Related Genes
After merging gene matrix of GTEx and TCGA datasets based on the Ensembl IDs, differentially expressed genes (DEGs) between pancreatic cancer and normal pancreas specimens were screened via the Linear Models for Microarray Data (limma) package in R (13). The screening criteria were as follows: |log fold change (FC)| > 2 and false discovery rate (FDR) < 0.05. Two thousand four hundred and ninety-eight immune-related genes were obtained from the IMMUPORT database (https:// www.immport.org/home). Following integration of DEGs and immune-related genes, abnormally expressed immune-related genes were identified for pancreatic cancer.

An Immune-Related Gene Signature Construction
For differentially expressed immune-related genes, univariate Cox regression analysis was conducted using TCGA dataset. Genes with p < 0.001 were considered related to pancreatic cancer prognosis. Then, candidate genes were screened via multivariate Cox regression analysis. Based on them, a risk score was established according to the following formula: risk score = β 1 x 1 + β 2 x 2 + . . . + β i x i (where βi indicates the coefficient of gene i, and x i indicates the expression level of gene i). The risk score of each patient was calculated and all patients were separated into high-and low-risk groups in accordance with the median value. Kaplan-Meier survival analysis was conducted between high-and low-risk groups through the survival package in R. The Receiver Operating Characteristic curves (ROCs) for overall survival were drawn utilizing the survivalROC package in R (14). Areas under the curves (AUCs) were calculated for detection of the efficacy to predict survival for the signature and other clinical features (age, gender, grade, pathologic T, pathologic N, and stage). Univariate and multivariate cox regression analyses were utilized to assess whether risk score could be independently predictive of patients' survival. Hazard ratio (HR), 95% confidence interval (CI) and p were calculated. HR > 1 indicated risk factors and HR < 1 indicated protective factors.

Functional Annotation Analysis
Gene ontology (GO) including biological processes (BP), molecular functions (MF) and cellular components (CC) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotation analyses of the survival-related genes were presented for exploring underlying biological functions via the clusterProfiler package in R (15). Terms with adjusted p < 0.05 were significantly enriched.

Protein-Protein Interaction (PPI) Analysis
The physical and functional interactions of proteins from the immune-related gene signature were analyzed via the Search Tool for the Retrieval of Interacting Genes Database (STRING) database (http://string-db.org/) (16). The degree of nodes was then calculated.

Gene Expression Profiling Interactive Analysis (GEPIA)
The expression of MET, OAS1, and OASL in pancreatic cancer and normal samples was confirmed using TCGA-GTEx projects via the online GEPIA database (http://gepia.cancer-pku.cn/ index.html) (17). Furthermore, overall survival (OS) and diseasefree survival (DFS) analyses were generated for high and low expression of MET, OAS1, or OASL groups among pancreatic cancer patients.

Immunohistochemistry and Immunofluorescence
Immunohistochemistry of MET, OAS1, and OASL proteins in pancreatic cancer and normal pancreas specimens was obtained from The Human Protein Atlas (https://www.proteinatlas. org/). Furthermore, their immunofluorescence images were also downloaded.

Immune Infiltration Analysis
The correlations between expression or copy number (armlevel deletion, diploid / normal, arm-level gain, and high amplification) of MET, OAS1, and OASL and the abundance of immune infiltrates composed of B cell, CD4 + T cell, CD8 + T cell, neutrophil, macrophage, and dendritic cell were assessed using the TIMER database (https://cistrome.shinyapps.io/timer/ (18)).

Statistical Analysis
All analyses were presented utilizing R version 3.5.2. Spearman correlation analysis was used to assess the associations between genes from the gene signature. The strengths of correlations were determined as follows: 0-0.39: weak; 0.40-0.59: moderate; 0.60-0.79: strong; 0.80-1.0: very strong. The differences in gene expression between two subgroups were calculated using the Wilcox test. P < 0.05 indicated statistical significance.

Construction of a Prognostic
Immune-Related Gene Signature for Pancreatic Cancer From TCGA-GTEx datasets, 1,737 DEGs with |log FC| > 2 and FDR < 0.05 were identified for pancreatic cancer (n = 178) than normal pancreas specimens (n = 165), which were composed of 962 down-and 775 up-regulated genes listed in Supplementary Table 1. Among all DEGs, 229 were immune-related genes ( Figure 1A; Supplementary Table 2). In Figure 1B, 50 abnormally expressed immune-related genes were significantly related to prognosis of pancreatic cancer. Using multivariate Cox regression analysis, 21 candidate genes were identified for constructing a prognostic immune-related gene signature ( Table 2). The risk score was calculated for each patient by combining coefficient and expression level. One hundred and seventy-eight patients were separated into high-and lowrisk groups in accordance with the median value. Those in the high-risk group exhibited an unfavorable prognosis (p = 1.499e-14; Figure 1C). To assess whether the risk score accurately and sensitively predicted patients' survival, ROCs were established. In Figure 1D, the AUC of the risk score for overall survival was 0.833, which was much higher than other clinicopathological factors, suggesting that the risk score possessed high accuracy in predicting survival.
The Immune-Related Gene Signature as an Independent Prognostic Factor for Pancreatic Cancer One hundred and seventy-eight pancreatic cancer patients were ranked according to their risk scores (Figure 2A). As the risk scores elevated, the number of dead patients gradually increased ( Figure 2B). Figure 2C depicted the expression patterns of these 21 genes in high and low-risk groups.  Figure 2D). Multivariate cox regression analysis revealed that risk score was an independent predictive factor for prognosis of pancreatic cancer ( Figure 2E).

Enrichment Analysis for Survival-Related Genes
We further probed biological functions of the survival-related genes. In Figure 3A, these genes were distinctly associated with regulation of migration of multiple cells such as epithelial and endothelial cells. Also, they could be involved in key cellular components like endoplasmic reticulum lumen and specific granule and possess different molecular functions such as receptor activity, cytokine binding and growth factor activity. As shown in KEGG enrichment analysis, these genes could participate in ErbB and proteoglycans in cancer pathways ( Figure 3B).

Interactions Between Genes From the Immune-Related Gene Signature
In Figure 4A, genes from immune-related gene signature were all abnormally expressed in pancreatic cancer than normal samples (all p < 0.001). Except for IL22RA1 and PAK3, most of them were up-regulated in pancreatic cancer. We further calculated their correlations at the expression levels using Spearson analysis, as shown in Figure 4B. S100A16 was strongly correlated to PAK3 (r = −0.68), SDC4 (r = 0.64), PLAUR (r = 0.74), MET (r = 0.66), and TMSB10 (r = 0.6). PAK3 had strong correlations with PLAUR (r = −0.64) and TMSB10 (r = −0.61). SDC4 exhibited a strong association with MET (r = 0.7). OASL was very strongly associated with OAS1 (r = 0.81). PLAUR showed a strong relationship with TMSB10 (r = 0.69). The PPI network confirmed the closely interactions between them ( Figure 4C). Figure 4D listed the   degree of each node in the network. We found that SSP1 had the highest degree.

Up-Regulation of MET, OAS1, and OASL Is Associated With Poor Clinical Outcomes of Pancreatic Cancer
Among 21 genes from the immune-related signature, MET, OAS1, and OASL were significantly associated with prognosis of pancreatic cancer patients. Patients with high MET expression were indicative of shorter DFS (p = 0.00044; Figure 5A) and OS (p = 0.00023; Figure 5B) time than those with low expression. Furthermore, high OAS1 (p = 0.047; Figure 5C) and OASL (p = 0.0072; Figure 5D) expression was distinctly related to poorer OS. MET (Figure 5E), OAS1 (Figure 5F), and OASL ( Figure 5G) were all up-regulated at the mRNA and protein levels. Immunofluorescence results demonstrated that MET (Figure 5H), OAS1 (Figure 5I), and OASL ( Figure 5J) were mainly distributed in cytoplasm and nucleus. This study further assessed whether MET, OAS1, and OASL expression was in association with clinical features. The data showed that MET expression was significantly higher in G3-G4 (p = 0.005; Figure 6A) and T3-T4 (p = 0.012; Figure 6B). Furthermore, higher OAS1 expression was detected in N1 (p = 0.019; Figure 6C) and T3-T4 (p = 0.006; Figure 6D). There was higher OASL expression in >65 (p = 0.048; Figure 6E) or T3-T4 (p = 0.009; Figure 6F) patients. These findings indicated that MET, OAS1, and OASL might be related to pancreatic cancer progression.

DISCUSSION
This study constructed an immune-related gene signature for predicting clinical outcomes of pancreatic cancer patients. We identified three key genes (MET, OAS1, and OASL) that were all up-regulated in pancreatic cancer and indicated an unfavorable prognosis. Following multiple dataset verification, the three genes might be promising therapeutic targets, which were worthy of further research. Despite the TNM stage system as an efficient tool to detect tumor stage, there is discrepancy in prognosis based on TNM stage (19). Thus, the efficacy of TNM stage is limited. Gene-based markers have been widely explored for pancreatic cancer in recent years (20). Recently, several prognosis-related gene signatures have been established for pancreatic cancer (21)(22)(23). For example, Zhuang et al. developed a prognosisrelated lncRNA signature for pancreatic cancer (24). Following comparing other clinical risk factors, the signature exhibited better predictability. Wu et al. constructed 9-gene signature for prediction of survival time of pancreatic cancer patients (25). There is still a lack of immune-related prognostic models. Herein, the immune-related 21-gene signature could accurately and sensitively predict survival time of pancreatic cancer patients. It performed better than other clinicopathological characteristics like age, gender, grade, N, T, and stage. This signature could be independently predictive of patients' prognosis.
The molecular mechanisms of highly aggressive behaviors remain unknown. We analyzed biological functions of survivalrelated genes. These genes could regulate migration of multiple cells such as epithelial and endothelial cells and were involved in key cellular components like endoplasmic reticulum lumen and specific granule and possess different molecular functions such as receptor activity, cytokine binding, and growth factor activity. These data were indicative that these genes were involved in tumor progression. For example, CCN1/Cyr61 secreted by pancreatic cancer cells may promote migration of endothelial cells (26). Nerve growth factors regulate CD133 functions, thereby promoting migration of pancreatic cancer cells (27). Targeting IL-1 and its receptor can prolong survival time in pancreatic cancer (28). We also found that these genes could participate in ErbB and proteoglycans in cancer pathways. It has been confirmed that dysregulated ErbB signaling promotes tumorigenesis for pancreatic cancer (29). Hence, these survival-related genes might participate in pancreatic tumor progression.
Our study revealed that there were physical and functional interactions between 21 genes from the immune-related gene signature. Among them, MET, OAS1, and OASL were verified to be markedly up-regulated in pancreatic cancer and associated with poor clinical outcomes of patients. Furthermore, MET expression was significantly correlated to infiltration of B cell, CD8 + T cell, CD4 + T cell, neutrophil, and dendritic cell. MET gene is located on chromosome 7q21-31. Changes in MET functions have been a hallmark of multiple cancers including pancreatic cancer (30). MET overexpression induces pancreatic cancer progression (31). Consistently, MET up-regulation was markedly correlated to tumor grade and T stage (32). Dysregulated MET functions correlate with aggressive phenotypes. Consistent with previous research, MET is involved in the crosstalk between tumor cells and tumor microenvironment (30). Thus, targeting MET has been considered as an adjuvant therapy in pancreatic cancer. OAS1 and OASL, 2 ′ -5 ′ -oligoadenylate synthetases, are interferon-induced antiviral enzymes. We found that OAS1 expression exhibited significant associations with neutrophil and dendritic cell infiltration. OASL expression showed a significant correlation with neutrophil infiltration. Consistently, Zhang et al. also found that OAS1 and OASL were correlated to neutrophil cell infiltration in breast cancer (33). Thus, MET, OAS1, and OASL were distinctly correlated to tumor immune microenvironment.
Collectively, we identified three immune-related prognostic genes MET, OAS1, and OASL, which could be promising therapeutic targets as well as prognostic markers for pancreatic cancer. However, several limitations of this study should be pointed out. First, the biological functions of MET, OAS1, and OASL such as their interactions with immune cells in tumor microenvironment need to be explored in vitro experiments. Second, their prognostic values should be verified in prospective research.

CONCLUSION
In this study, we established an immune-related 21-gene signature for prediction of pancreatic cancer prognosis. This signature could be accurately and independently predictive of  patients' survival. Among these genes, MET, OAS1, and OASL were validated to be up-regulated in pancreatic cancer and associated with unfavorable prognosis of patients. Also, there were closely interactions between them and immune infiltrates. Thus, MET, OAS1, and OASL could be potential therapeutic markers in pancreatic cancer.

AUTHOR CONTRIBUTIONS
HL conceived and designed the study. CZ, FN, and PH conducted most of the experiments, data analysis, and wrote the manuscript. YZo, YZh, YL, and HF participated in collecting data and helped to draft the manuscript. All authors reviewed and approved the manuscript.