Identification of hub genes and their novel diagnostic and prognostic significance in pancreatic adenocarcinoma

Objective: The main reasons for the poor prognoses of pancreatic adenocarcinoma (PA) patients are rapid early-stage progression, advanced stage metastasis, and chemotherapy resistance. Identification of novel diagnostic and prognostic biomarkers of PA is therefore urgently needed. Methods: Three mRNA microarray datasets were obtained from the Gene Expression Omnibus database to select differentially expressed genes (DEGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses for hub genes were performed using DAVID. Correlations between expression levels of hub genes and cancer-infiltrating immune cells were investigated by TIMER. Cox proportional hazard regression analyses were also performed. Serum hub genes were screened using the HPA platform and verified for diagnostic value using ELISAs. Results: We identified 59 hub genes among 752 DEGs. GO analysis indicated that these 59 hub genes were mainly involved in the defense response to viruses and the type I interferon signaling pathway. We also discovered that RSAD2 and SMC4 were associated with immune cell infiltration in the PA microenvironment. Additionally, DLGAP5 mRNA might be used as an independent risk factor for the prognoses of PA patients. Furthermore, the protein encoded by ISG15, which exists in peripheral blood, was validated as a potential diagnostic biomarker that distinguished PA patients from healthy controls (area under the curve: 0.902, 95% confidence interval: 0.819–0.961). Conclusions: Our study suggested that RSAD2 and SMC4 were associated with immune cell infiltration in the PA microenvironment, while DLGAP5 mRNA expression might be an independent risk factor for the survival prognoses of PA patients. Moreover, ELISAs indicated that serum ISG15 could be a potential novel diagnostic biomarker for PA.


Introduction
Pancreatic adenocarcinoma (PA) is a highly lethal tumor with a poor prognosis. Despite recent advances in diagnosis and treatment, its 5-year survival rate remains at 8% 1 . PA is predicted to be the second leading cause of cancer-related mortality by 2030 2 . Many PA patients postoperatively develop recurrences typically within 2 years and die within 5 years of recurrence [3][4][5] . Therefore, new requirements for diagnosis and treatment strategies have been proposed to improve the survival time of PA patients.
Recently, the increased use of immunotherapies has helped in the treatment of advanced PA patients. However, with feedback from clinical trials, PA is considered to be a type of "immune cold" tumor 6 . The outcome of many PA patients has not been significantly improved by immunotherapy, which poses a new challenge for immunotherapy strategies 7 .
Oncolytic viruses (OVs) are emerging immunotherapeutics for many cancers. Tumor-targeted OVs are usually engineered viruses that can selectively infect, replicate in, and lyse tumor cells, and also induce immune responses specifically against tumor cells 8 . OVs can also change the immunogenicity of the tumor microenvironment from an immunosuppressed (immunologically "cold") state to an immunoactivated (immunologically "hot") state 9 . The oncolytic herpesvirus, talimogene laherparepvec, was approved by the United States Food and Drug Administration and European Medicines Agency for advanced melanoma treatment in 2015. In addition, several OVs have been used in clinical trials of PA 10,11 . However, studies have found that some PA patients develop drug resistance due to the type I interferon (IFN) responses of the host that restrict viral replication 12 . The products of IFN-stimulated genes (ISGs) have been reported to contribute to the resistance of PA cells to OVs [13][14][15] . There is therefore an urgent need to identify specific defects in genes involved in IFN signaling pathways as potential biomarkers to classify candidate patients for effective oncolytic therapy.
The molecular study of carcinogenesis and progression is very important for patient therapy, diagnosis, and prognoses 16,17 . Using studies of the carcinogenesis and progression of PA, increasing attention has been directed toward abnormal expressions and mutations of genes [18][19][20][21] . In recent years, many studies have screened cancer-related differentially expressed genes (DEGs) by microarray technology. Furthermore, a microarray technique has been used to detect PA cells infected with an OV, and it was found that the expression levels of certain genes were significantly altered before and after viral infection, so these genes could be used as potential biomarkers to predict the effectiveness of oncolytic therapy 22 . However, the clinical sample size for PA is smaller than that of other common cancers 23 , leading to controversial results between different microarray chip analyses 24 . Thus, the integration of microarray technology and bioinformatics analysis technology is useful in identifying novel DEGs and functional pathways related to the occurrence and development of PA, as well as clinical responses to therapies [25][26][27] .
In the current study, we identified 752 DEGs between PA tissues and noncancerous tissues extracted from the Gene Expression Omnibus (GEO) database, among which 59 DEGs were identified as hub genes using Protein-Protein Interaction (PPI) network analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were used to characterize the hub genes at the molecular level. These hub genes were mainly enriched in the defense response to virus and the type I IFN signaling pathways. Furthermore, the correlation between hub genes and immune cell infiltration in the microenvironment was investigated by the Tumor Immune Estimation Resource (TIMER) webserver to provide new insight for OV therapy. More importantly, we screened potential biomarkers related to the diagnosis and prognosis of PA from 59 hub genes.

Microarray data
The mRNA microarray datasets, GSE15471, GSE16515, and GSE71989, based on the GPL570 platform, (HG-U133_Plus_2 Affymetrix Human Genome U133 Plus 2.0 Array) (Affymetrix, Santa Clara, CA, USA) were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo), which is a public functional genomics data repository containing microarray-and chip-based data of gene profiles 28

Identification of DEGs
The online GEO2R tool (http://www.ncbi.nlm.nih.gov/geo/ geo2r/) was used to screen DEGs between PA and normal pancreatic tissues. During analysis using GEO2R, the Benjamini & Hochberg (false discovery rate) method was used to adjust the P-values to balance the results of statistically significant genes and the limitations of false positives. Furthermore, a |log fold change (FC)| > 1 and adjusted P-value < 0.01 were further set as the cut-off criteria for screening DEGs.

GO functional and KEGG pathway enrichment analyses of DEGs
GO functional analysis is a useful bioinformatics tool for annotating genes and analyzing characteristic biological attributes for high-throughput genome or transcriptome data 32 . GO terms contain 3 main categories: biological process (BP), cellular component (CC), and molecular function (MF). KEGG (http://www.kegg.jp), an online database resource, including defined and associated gene sets and their pathways, was used for identifying high level functions and biological systems 33 . The GO functional analyses and KEGG pathway enrichment analyses of DEGs were generated by using Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8, http://david.ncifcrf.gov/) 34 . P < 0.05 and false discovery rate ≤ 0.01 were set as the filter criteria.

Visualized PPI network construction and hub gene selection
Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, http://string-db.org) 35 was used to predict a PPI network and evaluate the comprehensive information of proteins. In this study, a PPI network of all candidate DEGs was generated using the STRING website, and each interaction between proteins with a reliability threshold of a combined score of > 0.4 was considered statistically significant. The PPI network was further visualized by Cytoscape (version 3.7.1, http://www. cytoscape.org/) 36 . The Molecular Complex Detection (MCODE, version 1.5.1) plugin of Cytoscape was used to obtain the most densely connected module (the most significant clustering module) in the visible PPI network. The filter criteria for the most densely connected module were set as follows: "MCODE scores > 5", "node score cut-off = 0.2", "k-score = 2", "degree cutoff = 2", and "Max depth = 100". The hub genes were selected with degrees ≥ 20 in the most densely connected module.

TIMER database analysis
TIMER is a database that provides a deconvolution algorithm for systematic analysis to investigate the correlation between gene expressions and infiltrating immune cells, including B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells 37 . Correlations between the expression levels of hub genes and the infiltrating immune cells were analyzed using correlation modules. Moreover, the association between hub genes and tumor purity was also determined in the microenvironment of PA. The scatterplots were displayed by the TIMER web tool for correlation visualization. P < 0.05 was considered to indicate a statistically significant correlation.

Prognostic analyses of hub genes
The correlations of hub gene expression with overall survival (OS) and disease-free survival (DFS) were displayed by Kaplan-Meier curves from The Cancer Genome Atlas (TCGA) database in the Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) platform 38 . The median expression of the hub gene was set as the cut-off for demarcating high and low expression groups. A log-rank P-value < 0.05 was considered statistically significant. In addition, using the cBio Cancer Genomics Portal (cBioPortal, http://cbioportal. org) platform to identify the relationships between the alterations of hub genes and survivals of PA patients 39 , Kaplan-Meier curves were generated from the PanCancer Atlas of TCGA database. The following settings of genomic profiles were used for the analysis: "mutations", "putative copy number alterations from GISTIC", and "mRNA expression z-scores relative to diploid samples (RNA-seq V2 RSEM)". Complete samples (168 PA patients/cases) were selected for survival analyses, which included OS, DFS, progression-free survival (PFS), and disease-specific survival (DSS) analyses. Next, using univariate Cox proportional hazard analyses, risk factors were considered to be statistically significant with a P-value < 0.05, which was further analyzed by multivariate analysis to avoid confounding effects between each risk factor in the TIMER database 37 .

Diagnostic analysis of hub genes
To verify the mRNA expression level of the selected hub genes in pancreatic tissues, related published online data were obtained from TCGA and Genotype-Tissue Expression (GTEx) project datasets and analyzed using the GEPIA platform 38 . The mRNA microarray datasets were downloaded from the online Oncomine database (https://www.oncomine. org) to evaluate the diagnostic ability of hub gene expressions in tissues for PA. Circulating proteins encoded by these hub genes in blood were screened by the Human Protein Atlas (HPA) platform (https://www.proteinatlas.org/) and further verified for diagnostic value using ELISAs. To determine the diagnostic efficacies of hub genes for PA, MedCalc (version 18.2.1) was used to generate a receiver operating characteristic (ROC) curve and calculate the sensitivity and specificity of the area under the ROC curve (AUC). Hub genes that met the criteria (AUC > 0.5 and P-value < 0.05) were identified as potential diagnostic biomarkers.

Tissue microarray (TMA) and immunohistochemical analysis of the IFN-stimulated gene 15 (ISG15) protein
The TMA contains well-documented clinicopathological information, all of which was obtained with patient informed consent. The use of clinical samples was approved by the local Ethics Committee. TMA and ethical approval documentation were provided by Shanghai Outdo Biotech (Shanghai, China). The expression of ISG15 was analyzed using a TMA, which contained 35 cases of primary pancreatic ductal adenocarcinoma (PDAC) tissues and 22 cases of matched adjacent noncancerous tissues with 1 point for each tissue. Standard immunohistochemistry (IHC) procedures were performed on the TMA slide using ISG15 antibody (#sc-166755, 1:100; Santa Cruz Biotechnology). The staining intensity scores were no staining (score = 0), light brown (score = 1), brown (score = 2), and dark brown (score = 3), while the percentage scores of positive cells were determined as < 10% (score = 0), 10%-25% (score = 1), 26%-50% (score = 2), 51%-75% (score = 3), and > 75% (score = 4) 40 . The IHC ISG15 score was determined by multiplying the intensity score and the percentage score. IHC results were independently evaluated by 2 pathologists blinded to the clinical information. IHC ISG15 scores ≥ 6 in cancer tissues were defined as a "high expression level of ISG15."

TMA and immunohistochemical analyses of disks large-associated protein 5 (DLGAP5)
The details are presented in the supplementary materials and methods.

Validation of the diagnostic value of serum ISG15
PA patients were enrolled at the Department of Pancreatic Cancer of Tianjin Medical University Cancer Institute and Hospital from December 2019 to May 2020. Healthy controls were recruited at the Cancer Prevention Center of Tianjin Medical University Cancer Institute and Hospital. Participants with no evidence of pancreatic disease or other abnormalities based on the results of wellness check-ups including laboratory and radiological findings were enrolled as healthy controls. This study was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (bc2019104), and written informed consent was obtained from all PA patients and healthy controls.
A total of 78 serum samples were obtained from 40 PA patients and 38 healthy controls. The inclusion criteria were as follows: patients with histologically confirmed primary PA, PA patients admitted to the hospital for the first time without treatment, and PA patients and healthy controls aged > 18 years who did not have common viral infections. The exclusion criteria were as follows: PA patients or healthy controls taking immunotherapeutic drugs or antiviral drugs, PA patients or healthy controls with immune diseases, patients with both PA and other cancers, and healthy controls with a history of pancreatic disease.
The peripheral blood samples were centrifuged at 4 °C for 15 min at 3,000 rpm, and supernatant serum was stored at −80 °C until analysis. Serum ISG15 levels were analyzed according to the manufacturer's instructions using an ELISA kit (CSB-E12075h; Cusabio, Houston, TX, USA). The details are presented in the supplementary materials and methods.

Statistical analysis
Statistical analyses were performed using SPSS statistical software for Windows, version 25.0 (IBM, Armonk, NY, USA) and Prism 8.0 (GraphPad Software, San Diego, CA, USA). All quantitative data are expressed as the means and standard deviations from 2 or 3 separate experiments. For western blot analysis, statistical comparison between the 2 groups (PA cell vs. HPDE cell groups) was analyzed using Student's t-tests. For IHC analysis, categorical variables were analyzed using the chi-square test. The Kaplan-Meier method and log-rank tests were used to estimate the survival times of PA patients. For ELISAs, differences between the PA group and the healthy control group were assessed using the Mann-Whitney U test (continuous variables and nonparametric analyses). ROC curves were established to assess sensitivity, specificity, and AUCs with 95% confidence intervals (CIs). All P-values were 2-sided, and P < 0.05 was considered statistically significant.

Identification of DEGs in PA
To illustrate our study protocol, a flow chart is presented in Figure 1. The DEGs (1,815 in GSE15471, 1,775 in GSE16515, and 4,007 in GSE71989) were extracted from standardized microarray analyses. Three DEG datasets were then uploaded to a Venn diagram web-based tool (http://bioinformatics.psb. ugent.be/webtools/Venn/) to obtain the intersection set, which included 752 DEGs between PA tissues and noncancerous tissues (Supplementary Figure S1), with 621 upregulated genes and 131 downregulated genes.

GO functional and KEGG pathway enrichment analyses of 752 DEGs
As shown in Figure 2 and Supplementary Table S1, GO analysis revealed that 752 DEGs were significantly enriched in BPs, mainly including extracellular matrix (ECM) organization, type I IFN signaling pathway, cell adhesion, and collagen catabolic process. For MF, these DEGs were significantly enriched in integrin binding, collagen binding, ECM structural constituent, and laminin binding. For CC, the extracellular exosome, extracellular space, ECM, and proteinaceous ECM were significantly enriched. KEGG pathway analysis showed that these DEGs were significantly enriched in ECM-receptor interaction, amoebiasis, focal adhesion, and phagosome.

Visualized PPI network construction and module analysis
The visualized PPI network of 752 DEGs was generated and displayed using Cytoscape (Supplementary Figure S2). The most densely connected module was extracted from this PPI network (Figure 3 and Supplementary Figure S2). A total of 59 genes in this module were identified as hub genes with degrees ≥ 20. Geminin (GMNN) was downregulated, and other hub genes were upregulated. These genes were mainly enriched in the defense response to virus, response to virus, mitotic sister chromatid segregation, cell division, mitotic nuclear division, and type I IFN signaling pathways (Figure 4 and Supplementary Table S2). Two significantly enriched pathways in KEGG were herpes simplex infection and influenza A (Figure 4 and Supplementary Table S2).

The Hub genes related to immune cell infiltration were selected
We speculated that the results of GO and KEGG enrichment analyses of these 59 hub genes were related to the regulation of the immune system. It was therefore necessary to further understand the association of hub genes and immune cell infiltration in the PA microenvironment. Scatterplots and fitting curves of the correlations between hub genes and immune cell infiltration were generated by TIMER ( Figure 5 and Supplementary Table S3). The results suggested that the expressions of cytoskeleton associated protein 2 (CKAP2), guanylate binding protein 1 (GBP1), IFN gamma inducible protein 16 (IFI16), IFN regulatory factor 9 (IRF9), 2′-5′-oligoadenylate synthetase 2 (OAS2), radical S-adenosyl methionine domain containing 2 (RSAD2), structural maintenance of chromosome 4 (SMC4), and ubiquitin conjugating enzyme E2 L6 (UBE2L6) had significant correlations with the infiltration levels of all 6 immune cells (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells) (P-value < 0.05). Other hub genes were correlated with some of the above infiltrating immune cells. Additionally, RSAD2 was more closely associated with macrophages, neutrophils, and dendritic cells than with other infiltrating immune cells (partial correlation > 0. 3

Diagnostic capability evaluation of hub genes
To assess the mRNA expressions of 59 individual hub genes in PA and normal pancreatic tissues, we analyzed 8 datasets from Oncomine and found that the differential mRNA expression levels of the major histocompatibility complex, class I, A (HLA-A), IFN alpha inducible protein 27 (IFI27), ISG15, and karyopherin subunit alpha 2 (KPNA2) approximately correlated with the top 1% of genes in gene rank in the 3 datasets (Supplementary Figure S6). High mRNA expression levels of HLA-A, IFI27, ISG15, and KPNA2 in PA tissues vs. normal pancreatic tissues were considered statistically significant (P < 0.05) in 5 or 6 datasets (Supplementary Table S4). Next, 350 samples from TCGA and GTEx databases were used to analyze the mRNA expression levels of ISG15 in PA vs. normal tissues. The overexpression of ISG15 in PA tissues was significantly higher than that in normal tissues (P < 0.05) (Figure 6A)  dataset (P = 0.0063) 41 ; AUC = 1.000 with 100% sensitivity and 100% specificity in the Iacobuzio-Donahue dataset (P < 0.0001) 42 ; AUC = 0.920 with 80.0% sensitivity and 100% specificity in the Logsdon dataset (P < 0.0001) 43 ; AUC = 0.911 with 83.3% sensitivity and 100% specificity in the Pei dataset (P < 0.0001) 30 ; and AUC = 0.773 with 72.7% sensitivity and 100% specificity in the Segara dataset (P = 0.0302) 44 ] ( Figure  6B and 6C). Similarly, we evaluated the diagnostic value of HLA-A, IFI27, and KPNA2 mRNA expressions in the same way (Supplementary Figure S7, Figure S8, and Figure S9).

External validation of the protein levels of the hub genes
We then measured the protein levels of the hub genes that may have potential clinical value. We detected the protein levels    The rank for a gene is the median rank for that gene across each of the analyses. The P-value for a gene is its P-value for the median-ranked analyses.  Figure S11A). The protein levels of DLGAP5, KPNA2, and IFI27 in 2 kinds of PA cells (BxPC-3 and MIA PaCa-2) were higher than those in normal pancreatic cells (Supplementary Figure S11B). The protein levels of HLA-A and ISG15 in 3 kinds of PA cells were higher than those in normal pancreatic cells (Supplementary Figure S11B). Using the HPA platform to determine protein levels for these genes in clinical samples, we observed that only the protein encoded by ISG15 was expressed in tissues and secreted into the blood. We further detected the protein level of ISG15 in the TMA containing PDAC and normal adjacent pancreatic tissues ( Figure 7A) and found that staining for ISG15 was strong in PDAC tumor cells (Figure 7A, top), whereas in noncancerous tissues, the staining was weak or undetectable, except in islet cells that showed positive staining (Figure 7A, bottom). The levels of ISG15 protein in cancerous tissues were significantly higher than those in adjacent noncancerous tissues (Figure 7A, 7C,  7D). In all PDAC tissues, the number of cases with high levels of ISG15 with IHC scores greater than or equal to 6 accounted for 45.7% (Figure 7C, left), the number of cases with scores between 3 and 5 accounted for 17.2% (Figure 7C, left), and the number of cases with scores between 0 and 2 accounted for 37.1% (Figure 7C, left). In all adjacent noncancerous tissues, 100.0% exhibited low levels (IHC scores less than 6, Figure 7D), in which the number of cases with IHC scores between 0 and 2 accounted for 90.9% (Figure 7C, right), and those between 3 and 5 accounted for 9.1% (Figure 7C, right). However, the ISG15 protein level showed no significant correlation with other clinicopathological factors including age, gender, pathological grade, T (tumor) stage, N (node) stage, and M (metastasis) stage (Supplementary Table S5). The ISG15 protein levels were evaluated to differentiate cancerous tissues from adjacent noncancerous tissues using ROC curve analysis (AUC: 0.729, 95% CI: 0.594-0.838, sensitivity: 45.7%, specificity: 100.0%) ( Figure 7E). Next, we validated the protein levels of serum ISG15 in PA patients and healthy controls using ELISAs. Thirty-eight healthy participants were included as healthy controls (16 males and 22 females), and 40 PA patients (16 males and 24 females) were included in this study. The age range of PA patients was 36 years to 75 years, and the age range of healthy controls was 29 years to 57 years. Significant differences in serum ISG15 levels were found between the 2 groups (P < 0.001) ( Figure 7F). ROC curve analysis was performed to differentiate PA from healthy controls (AUC: 0.907, 95% CI: 0.819-0.961, sensitivity: 75.0%, specificity: 100.0%) ( Figure 7G). The results showed that the ISG15 protein could be a potential diagnostic biomarker of PA in both serum and tissues.

Discussion
In our study, the majority of 752 DEGs were enriched in the ECM-receptor interaction pathway. It is well-known that the ECM is one of the unique characteristics of PA. Anomalous proliferation and deposition of the ECM promote the formation of the microenvironment in PA, thus facilitating tumor progression, metastasis, and therapeutic resistance 45,46 . In addition, we found that the most densely connected module from the PPI network consisted of 59 hub genes with degrees ≥ 20.  follows: samples with IHC scores < 6 were defined as low expression, and scores ≥ 6 were defined as high expression. The percentage of cells staining positive was scored as follows: < 10% (score = 0), 10%-25% (score = 1), 26%-50% (score = 2), 51%-75% (score = 3), and 76%-100% sensitivity of the cells to OV infection 13,49,50 . These genes may be used as potential biomarkers to predict the success of OV infection and replication 13 .
In the present study, we determined the associations of hub genes and immune cell infiltration in the microenvironment of PA using the online TIMER platform. We found for the first time that 8 hub genes (CKAP2, GBP1, IFI16, IRF9, OAS2, RSAD2, SMC4, and UBE2L6) were related to immune cell infiltration. GBP1 51 , IFI16 52 , IRF9 51 , OAS2 47 , and RSAD2 53 were demonstrated to be restriction factors in some OV infections or replications. Because these genes were mainly enriched in the defense response to viruses (Supplementary Table S2), we speculated that high expression levels of these genes may have led to resistance of PA cells to OV treatments, which requires confirmation in future studies. RSAD2 is an ISG that is a broadly acting effector of the type I/II IFN antiviral response against several enveloped viruses, such as herpes virus, influenza A virus, and measles virus (MV) [53][54][55] . Kurokawa et al. 53 reported that RSAD2 was a restriction factor of the oncolytic MV that inhibited the release of MV in infected ovarian cancer cells (SR-B2). Human glioblastoma cells were infected by attenuated MV as a kind of oncolytic virotherapy, which activated RSAD2 expression in tumor cells 47 . They suggested that RSAD2 activation could be a predictive biomarker for the response to oncolytic MV therapy. However, the role of RSAD2 in OV therapy for PA remains unclear and should be further investigated in future studies. The structural maintenance of chromosome 4 protein is encoded by SMC4, which is a structural maintenance protein of the chromosome (SMC) family member. High expression of SMC4 is related to tumor size, dedifferentiation, advanced stages, and vascular invasion of liver cancer 56 . In glioma and colorectal cancer, upregulated expression of SMC4 promoted the tumor cell growth rate, migration and invasion [57][58][59] . Overexpression of SMC4 is an independent prognostic factor in lung adenocarcinoma 60 . High expression of SMC4 predicted worse prognosis of pancreatic ductal adenocarcinoma patients 61 .
In the present study, we first observed that SMC4 was associated with immune cell infiltration in the microenvironment of PA. Moreover, SMC4 was correlated with the OS and DFS of PA patients. Nevertheless, the role of SMC4 in the immune microenvironment of PA remains unclear and should be further investigated in future studies.
Our results indicated that alterations in ANLN, ASPM, CDK1, CEP55, DLGAP5, and NUSAP1 were significantly related to the survival of PA patients. DLGAP5 is an important mitotic spindle protein, and the upregulation of DLGAP5 contributes to the proliferation, migration, and invasion of multiple cancers [62][63][64] . Overexpression of the DLGAP5 protein or increased mRNA expression are associated with an unfavorable prognosis in multiple cancers, such as colorectal cancer 64 , ovarian cancer 65 , and lung cancer 66,67 . However, until recently, only 1 study has been performed to identify DLGAP5 mRNA expression as a potential biomarker for the diagnosis and prognosis of PA by bioinformatics analysis 68 . To the best of our knowledge, there are no reports of detecting DLGAP5 protein levels in PA tissues. In the present study, we found that alterations in DLGAP5 were related to worse OS, DFS, PFS and DSS in PA patients. More importantly, DLGAP5 mRNA expression might be an independent risk factor for the prognoses of PA patients, using bioinformatics analyses. To determine whether the protein level of DLGAP5 has prognostic value for PA, we used immunohistochemical methods to detect the protein level of DLGAP5 in the TMA of PA patients with prognostic information. No positively stained cells were found in any adjacent noncancerous tissues (an IHC DLGAP5 score of 0 was defined as a negative result), while some cancerous tissues had a small number of positively stained cells (an IHC DLGAP5 score of 1 or 2 was defined as a weakly positive (score = 4). (B) Representative images of PDAC tissues with no staining/weak/moderate/strong ISG15 expression magnified 200-fold. Scale bars: 100 μm. The staining intensity was scored as follows: no staining (score = 0), weak (score = 1), moderate (score = 2), and strong (score = 3). (C) The proportion of cancerous (Figure 7C, left) and adjacent noncancerous (Figure 7C, right) tissues with ISG15 protein levels based on IHC scores (0 to 12). (D) Correlation between the protein level of ISG15 and histological type. IHC scores ≥ 6 were considered high levels of ISG15. A P-value < 0.05 was considered to be statistically significant. (E) An ROC curve of ISG15 protein levels [high levels of ISG15 (IHC scores ≥ 6) and low levels of ISG15 (IHC score < 6)] was used to differentiate cancerous tissues (n = 35) from adjacent noncancerous tissues (n = 22) (P < 0.001). (F) High protein levels of serum ISG15 in PA patients (n = 40) compared with those in healthy controls (n = 38) (P < 0.001). (G) A ROC curve of the serum ISG15 levels was used to differentiate PA patients from healthy controls. PDAC, pancreatic ductal adenocarcinoma; IHC, immunohistochemistry; TMA, tissue microarray; n, number, ROC, receiver operating characteristic; AUC, area under the ROC curve; CI, confidence interval. Figure S12A, S12B). Although we found that weakly positive expression of DLGAP5 protein was associated with higher pathological grade and shorter OS of PDAC patients in this study (P < 0.05) (Supplementary Table  S6 and Figure S12C), we believe that whether the protein level of DLAGP5 has prognostic value needs to be further confirmed by increasing the number of clinical specimens. We will collect more clinical samples in the future to verify the protein level of DLGAP5 in PA and normal adjacent pancreatic tissues and to explore its prognostic value. In addition, studies on the regulatory mechanism of DLGAP5 in PA have rarely been reported, which should be conducted and validated in future investigations.

result) (Supplementary
ISG15 encodes a 15-kDa ubiquitin-like protein, which is induced by IFN-α and IFN-β. KRAS-associated phenotypes of PA cells are reversed by knocking down the ISG15 pathway proteins (free ISG15 and ISG15 conjugates), which by itself or in combination with anti-PD-1 treatment may contribute to improved survival for patients with PA 69 . Another study reported that extracellular ISG15 maintained the stem celllike features of PA cells, such as clonogenicity, invasiveness, and spheroid formation, by an autocrine mode of action 70 . Using data analysis of the Oncomine database, we showed that ISG15 mRNA expressions in tissues may have potential value in the diagnosis of PA. In terms of the protein level of ISG15 in clinical samples, little is known about its diagnostic value in PA. In our study, the differential protein levels of ISG15 were measured between PA and adjacent noncancerous tissues by IHC. The protein level of serum ISG15 was further validated by ELISA, and the diagnostic role of serum ISG15 was identified by ROC analysis.

Conclusions
Using bioinformatics analyses, our study identified DLGAP5 mRNA expression as a potential independent risk factor for prognoses of PA patients, while showing that RSAD2 and SMC4 were related to immune cell infiltration. The protein levels of ISG15 in PA tissues were significantly higher than those in adjacent noncancerous tissues, which were validated by IHC. Additionally, the differential expression of serum ISG15 was validated by ELISAs to show that it could successfully distinguish PA patients from healthy controls by ROC analysis. However, there are several limitations in the current study. First, based on the present data from bioinformatics analyses, functional experiments are needed to verify these results. Moreover, we observed that serum ISG15 was a potential novel non-invasive biomarker for the diagnosis of PA, but this was based on a single cohort with a small sample size. Future studies using larger independent cohorts are needed to validate our findings, and the control group should involve chronic pancreatitis patients.