A Novel Immune-Related Prognostic Signature Predicting Survival in Patients with Pancreatic Adenocarcinoma

Pancreatic adenocarcinoma (PAAD) carries the lowest survival rate of all major organ cancers, which is of dismal prognosis and high mortality rate. Thus, the present study attempted to identify a few novel prognostic biomarkers and establish an immune-related prognostic signature which could predict the prognosis of PAAD. Four prognostic immune-related genes (IRGs) including S100A6, S100A10, S100A16, and SDC1 were screened by differentially expressed gene (DEG) identification and weighted gene coexpression network analysis (WGCNA). Subsequent analysis proved the high expression of these IRGs in PAAD tissues, suggested by TCGA-PAAD data, merged microarray-acquired dataset (MMD), GEPIA, and Oncomine webtool. By using MMD and TCGA-PAAD data, S100A6 (MMD: AUC = 0.897; TCGA: AUC = 0.843), S100A10 (MMD: AUC = 0.880; TCGA: AUC = 0.780), S100A16 (MMD: AUC = 0.878; TCGA: AUC = 0.838), and SDC1 (MMD: AUC = 0.885; TCGA: AUC = 0.812) exhibited excellent diagnostic efficiency for PAAD. By conducting connectivity map (CMap) analysis, we concluded that three molecule drugs (sulpiride, famotidine, and nalidixic acid) might have worked in the treatment of PAAD. Then, an immune-related prognostic index was constructed, which was validated as an independent prognostic factor for PAAD patients (P=0.004). We further constructed a nomogram by using this immune-related signature and age, the prognostic value of which was validated by using concordance index (C-index = 0.780) and area under curve (AUC = 0.909). Moreover, the immune-related prognostic signature was associated with response to anti-PD-1/L1 immunotherapy. To sum up, four IRGs were screened out and verified to be novel immune-related prognostic biomarkers in PAAD. Besides, sulpiride, famotidine, and nalidixic acid might be potential choices in the treatment of PAAD. An immune-related signature was established to show great potential for prognosis prediction for PAAD, independently, which might guide more effective immunotherapy strategies. A nomogram is further established by using this immune-related prognostic index, which might contribute to more effective prognosis prediction in PAAD patients.


Introduction
Pancreatic adenocarcinoma (PAAD) mainly refers to a malignant tumor derived from pancreatic duct epithelial cells and follicular cells, which is a common malignant tumor of digestive system [1]. According to the data of the National Comprehensive Cancer Network Center in 2020, about 540,000 people were diagnosed with pancreatic cancer in the United States in 2019, and the number of deaths due to pancreatic cancer was about 430,000 [2]. e incidence rate of pancreatic cancer is the top ten among the leading cancers in the United States, ranking third in the cancer-related death causes, as the American Cancer Society reported [2]. Pancreatic ductal adenocarcinoma (PDAC) is the most common type of pancreatic adenocarcinoma, which comes from pancreatic duct epithelial cells, accounting for 80%∼90% of all pancreatic cancer patients [3]. Some studies indicate that the mortality of PDAC is tightly related to incidence [4]. And what is worse is the mortality of PDAC ranks the fourth commonest reason of cancer-related death worldwide [5]. Nowadays, the fundamental treatment principle is still based on surgery, combining with comprehensive treat methods such as radiotherapy and chemotherapy. Despite great efforts, PDAC is still of poor survival, 5-year survival of which is less than seven percent [6,7]. Because of the poor prognosis of PAAD, there might be a great need of novel prognostic biomarker exploring.
As a kind of cancer therapy, immunotherapy could fight cancer based on immune system [8]. At present, there is a growing realization that immunotherapy might show potential to treat tumors effectively and safely [9][10][11]. Bioinformatics have made significant progress recent years, which promoted the increasing use of mining of public databases in cancer biomarkers identification. With theory researches going deep, immune-related genes (IRGs) have shown potential among development and immunotherapy of several cancers [12][13][14]. Chen et al. have developed a prognostic signature for head and neck squamous cell carcinoma via IRGs, which could effectively distinguish the prognosis, the molecular and immune characteristics, and immune benefit from immune checkpoint inhibitor therapy [15]. ese researches determined that immune-related predictive biomarkers could enhance immunotherapy efficacy [16,17]. Screening of prognostic biomarkers via IRGs has without question been immune therapy of cancer research hotspot. Lack of prognostic biomarkers related to the tumor immune microenvironment for PAAD patients stimulates us to explore some immune-related biomarkers in PAAD, which might guide appropriate therapy tips to improve the therapeutic efficacy in PAAD.
So far as we know, this research might be the first one to screen immune-related prognostic biomarkers for PAADs straight by bioinformatics. We firstly collected nine datasets from Gene Expression Omnibus (GEO) and e Cancer Genome Atlas (TCGA) database. en conducting WGCNA via immune-related genes in PAAD. Hub genes reaching the standards have been screened out from the key module. Furthermore, 22 overlapped DEGs were selected, four of which were common with hub genes. us, the four IRGs were considered to be initial prognostic factors and further validated using other methods. After that, a risk signature via the four IRGs was constructed, which might be a novel prognosis prediction tool for PAAD and guide more effective immunotherapy strategies. In addition, we developed a nomogram comprehensively considering the risk signature, clinical characters. e nomogram established for forecasting the survival rate of PAAD could give some beneficial guidance for clinical application.

Dataset and Immune-Related Gene Collection.
e research steps in this study are described in Figure 1, showing the identification and validation of immune-related prognostic biomarkers in PAAD. PAAD microarray data were firstly downloaded from TCGA database (https:// genomecancer.ucsc.edu/). e microarray matrix were displayed as count number. We removed tumor samples without complete clinical information from subsequent analysis. In other words, all the samples left contained completed survival information and essential clinical factors (age, gender, tumor grade, and tumor stage). TCGA-PAAD data standardization including normalization and log2 transformation were in process having the aid of R package "DEseq.2." 4 normal samples and 177 PAADs were used for subsequent process. en four independent GEO datasets (GSE15471 [18], GSE16515 [19], GSE22780, and GSE32676 [20]) were collected through GEO database (http://www.ncbi.nlm.nih. gov/geo/). Because these datasets used the same annotation platform (GPL570), we merged them to develop a merged microarray-acquired dataset (MMD). We firstly downloaded the raw data for the four datasets. Secondly, we normalized the raw data using RMA-normalization method, based on package "affy" [21] in the R software. en, we preprocessed, merged, and ComBat-adjusted the four datasets based on R package insilico-Merging. A PAADspecific, MMD was generated though the above steps. Finally, we used the GPL570 annotation files for annotation of probes. After finishing these procedures, 70 normal tissues and 108 PAAD tissues were collected and further used in this research.
Besides, with the aim of validation of immune-related prognostic biomarkers, three independent GEO datasets (GSE21501 [22], GSE28735 [23], and GSE71729 [24]) with complete survival information were retrieved from GEO database. Raw expression data of GSE28735 was firstly retrieved from GEO database and further normalized and transformed via "affy" package, as previous done. For GSE21501, we used locally weighted linear regression (Lowess) for normalization. For GSE71729, nonnegative normalization was conducted. en we merged the three datasets to develop a merged microarray-acquired survival dataset (MMSD) as we previous did. Totally 269 PAADs were included in MMSD. e detail information of the eight datasets were showed in Table 1.
IRGs were retrieved from ImmPort database (https:// www.immport.org). We obtained 2,499 IRGs from this database. e 1,650 overlapped genes with the TCGA-PAAD gene list were selected for next-step study.

WGCNA to Screen Key Module.
Firstly, the expression matrix of the 1,650 IRGs was checked via two approaches (goodSamplesGenes and sample-network-method). Outliers were further identified with this cut-off criterion of Z.Ku < −2.5 (Z.ku � (ku-mean(k))/(sqrt(var(k)))). We removed the substandard samples from further analysis. en a coexpression network was established via "WGCNA" [25]. Branch cutting methods was conducted to classify IRGs into gene modules [26]. Some important parameters in branch cutting methods were set including minClusterSize � 30, and deepSplit � 2. After splitting IRGs into gene modules, a cut line (correlation ≥ 0.75) for combining high related modules was selected, which was realized by measuring the dissimilarity of module eigengenes (MEs). With the aim of screening hub modules associated with disease status (PAAD or normal, the trait interested us most), we firstly calculated the Gene Significance (GS). Furthermore, module significance (MS) was calculated based on GS. MS was measured as the average GS of all the genes in a module. After finishing the above-mentioned processes, the most correlated module was selected, which might be a key   module. Furthermore, we regarded genes with the cut-off criterion (|cor.geneModuleMembership| >0.8 and |cor.ge-neTraitSignificance| >0.2) as key genes in WGCNA, which were included for next-step analysis.

Differentially Expressed Immune-Related Gene
Identification. Based on TCGA-PAAD data (1,650 IRGs), we firstly screened out differentially expressed genes (DEGs) between normal tissues and PAADs by R package "edgeR" [27]. Similarly, based on "limma" [28], DEGs were identified by using the expression matrix of 1,242 IRGs in MMD. Genes with adjusted P value <0.05 and |log2FC| ≥ 1.0 were regarded significantly different expressed. Finally, DEGs overlapping between DEGs from TCGA-PAAD and DEGs from MMD were chosen for further analysis.

Hub Gene Identification.
Combining WGCNA and DEG, hub gene were identified. Genes which could be found both in key biomarkers in WGCNA and DEGs were considered to be hub genes for further analysis. Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed via R package "clusterProfiler" [29] for functional annotation of hub genes. For the GO part, we just obtained the biological process (BP). We selected P < 0.05 as the standards to define significant BPs and KEGG pathway terms.

Small Molecule Drug Exploring.
We also attempted to explore some small molecule drugs which might be novel choices for PAAD therapy. Based on R package "edgeR", DEGs were firstly obtained based on the expression data profile of TCGA-PAAD data. DEGs were screened out using the same cut-off criterion we set before. Based on these DEGs, we conducted connectivity map (Cmap) (https:// portals.broadinstitute.org/cmap/) [30] analysis to explore small molecule drugs. A small molecule drug with P value <0.05 and |mean ≥ 0.50 was concluded to show strong potential in the treatment of PAAD.

Internal Validation of Hub Genes.
To explore and verify the prognostic value of key biomarkers, Gene Expression Profiling Interactive Analysis (GEPIA) webtool [31] was used to perform overall survival (OS) and disease-free survival (DFS) analyses (http://gepia.cancer-pku.cn/). Furthermore, the expression difference between normal samples and PAADs was also obtained via this webtool (to validate the result of DEG screening).

Verification of Expression of Hub Genes.
Besides exploration via GEPIA, the mRNA-level expression difference of hub IRGs was verified with the help of Oncomine database (https://www.oncomine.org/) [32]. Seven analysis were retrieved from Oncomine database and used for expression exploration. e mRNA-level expression of hub genes between normal tissues and PAAD tissues were also obtained via e Human Protein Atlas (HPA) database (https://www. proteinatlas.org/).

Prognostic Role of Hub Genes Exploration.
To validate the prognostic value of hub genes, survival analysis and receiver operating characteristic (ROC) analysis were conducted. Based on GSE21501 (n � 102), GSE28735 (n � 42), GSE57495 (n � 63), and GSE71729 (n � 125), PAADs were divided into high-and low-expression groups respectively according to the best calculated cut-off by R package "maxstat" from the beginning. OS analysis was next performed via R package "survival" [33]. Secondly, by using TCGA-PAAD and MMD large datasets, we plotted ROC curves to see if hub genes could distinguish PAADs and normal samples. By R package "pROC" [34], we did this analysis. Moreover, the area under curve (AUC) was also calculated. We considered a hub gene to have strong prognostic value when AUC > 0.75 in both the datasets.

Association between Hub Gene Expression and Immunocytes Exploring.
ere is a realization that immunocytes might act as independent predictors of survival in cancers. erefore, the association between IRGs and immunocytes were explored via TIMER (https://cistrome.shinyapps.io/ timer/) [35]. A hub gene with |correlation coefficient (cor) | ≥0.2 and P value <0.05 was thought to strongly relate to an infiltrating level of an immune cell type as previous did.

Functional Exploration of Hub Genes.
Gene set enrichment analysis (GSEA) might have researchers to understand the role of genes in biological behaviors. To explore the biological functions of hub IRGs, the median value of each hub gene was calculated via TCGA-PAAD data. 177 PAADs were divided into two groups immediately (highand low-expression groups). We selected "c2.cp.kegg.v7.4.symbols.gmt" as the annotated gene-set. After finishing these steps, we performed GSEA (http:// software.broadinstitute.org/gsea/index.jsp) [36], biological pathways of nominal P < 0.05, |ES| > 0.6, gene size (n) ≥ 20 and FDR <25% were considered significant.

Construction of an Immune-Related Prognostic
Signature. To comprehensively consider the prognostic value of hub genes, univariate Cox analysis of OS was firstly obtained to identify prognostic IRGs. P < 0.05 was set as the cut-off criteria to screen prognostic biomarkers. en relying on the expression levels and regression coefficient (Coef ) of prognostic biomarkers, we established an immune-related prognostic signature. Furthermore, we evaluated risk scores (RS) of PAADs as follows: Coef is defined as the regression coefficient of a prognostic biomarker and Exp betokens a prognostic biomarker expression level. We further evaluated the RSs of PAAD samples retrieved from TCGA-PAAD data and MMSD using the equation, in order to measure the prognostic value of the signature. PAADs were classified into two series (highand low-risk) in all datasets, respectively, according to the median RS in each dataset. Furthermore, we obtained OS analysis with the help of R package "survival". In addition, time-independent (1-, 3-, and 5-year) receiver operating characteristic (ROC) curves were also plotted by using R package "survivalROC" [37].

Cox Proportional Hazards Regression Analysis.
With the aim of prognostic value of hub genes validation, the risk score assessed by this immune-related signature and other essential clinical features (gender, age, tumor grade, and pathologic stage) from TCGA-PAAD data were selected for OS univariable Cox analysis. A factor of P value <0.05 was identified and further selected for multivariate Cox analysis. It will be sure if the immune-related signature was independent from the rest clinical factors for predicting OS of PAADs via this analysis. en package "forestplot" [38] in R software was used for visualization.

Establishing and Validating of a Nomogram.
With the aim of avoiding the overfitting problem, we conducted cross-validation before nomogram construction. en based on the immune-related prognostic index, a nomogram was constructed via R package "rms". Calibrate curves were drawn to test the nomogram, the 45°line in which was defined as the best prediction. In addition, we calculated the consistency index (C-index) between actual probability and predicted probability to further measure the prediction effectiveness of the nomogram. Via "pROC" in R software, we also conducted ROC analysis. Furthermore, time-dependent (1-, 3-, 5-year) ROC analysis was conducted. We immediately obtained decision curve analysis (DCA) via R package "rmda" [39] to explore the clinical application value of the nomogram. Survival analysis and GSEA were also performed to explore the survival difference and lucking function of the prognostic signature.
2.14. Immune-Related Prognostic Signature in the Role of Anti-PD-1/L1 Immunotherapy. Two immunotherapeutic cohorts including their clinical information were collected and included in this step. For the IMvigor210 cohort which contained advanced urothelial carcinoma, we retrieved the expression data displayed as count number via http://research-pub.Gene.com/imvigor210corebiologies. Relying on R package "DEseq2" [40], the count-numberbased matrix was transformed into TPM values. en totally 298 samples with intervention of atezolizumab (an anti-PD-L1 antibody) were collected and used for further exploration. Moreover, a dataset (GSE78220) [41] including metastatic melanoma samples treated with pembrolizumab (an anti-PD-1 antibody) was retrieved from GEO database. R package "limma" was used for the normalization, the FPKM value was also transformed into the TPM value. In total, 27 samples were collected and used for further exploration. e survival curves for prognostic analysis were generated by using R package "survival". And receiver operating characteristic (ROC) curves were also drawn and the area under the curve (AUC) were further calculated by using R package "pROC". Moreover, we also explored the association between the immune-related prognostic signature and immunocytes by using the Cell type Identification by Estimating Relative Subsets ff RNA Transcripts (CIBER-SORT) (https://cibersort.stanford.edu/). In addition, by using PACA-CA data with tumor status information (primary or metastatic) retrieved from ICGC data portal (https://dcc.icgc.org/), we also explored the relationship between the immune-related signature and tumor status.

Key Module Identification.
After removing 8 outliers, totally 173 PAADs were used for WGCNA ( Figure S1). In order to evaluate adjacencies, beta (β) � 9 (scale free R 2 � 0.84) was set as the soft-thresholding power ( Figure S2). IRGs were identified and further assigned to gene modules. In total, six modules were identified, which were showed in Figure 2(a). Genes with weak association with clinical trait were incorporated into the grey module and further removed in the present study. Furthermore, yellow module was chosen among the six modules because of the most positive correlation with disease status (P � 0.02, R 2 � 0.18, Figure 2(b)). As Figure 2(c) suggested, MM and GS of the yellow module (cor � 0.33, P � 0.023) showed significant relationship. We also found that the MS of the yellow module was the highest compared with the rest modules ( Figure 2(d)). erefore, we considered yellow modules as key module. Figure S3(a) showed the network heatmap based on these IRGs. In addition, the classical MDS plot ( Figure S3(b)) suggested that the six modules were independent from each other.  Table S1 showed the detail information of different expressed IRGs. We further conducted GO and KEGG pathway analysis for functional exploration. As shown in Figure 4(a), GO analysis indicated that IRGs were involved in negative regulation of endopeptidase activity, negative regulation of peptidase activity, negative regulation of proteolysis, regulation of endopeptidase activity, positive regulation of innate immune response, and regulation of peptidase activity. Furthermore, the selected IRGs were significantly associated with IL-17 signaling pathway, suggested by KEGG enrichment analysis ( Figure 4(b)).

Novel Choices for PAAD Treatment.
In order to provide some drugs to treat PAAD, we also performed CMap analysis. DEGs were firstly screened out by using MMD. Totally, 844 DEGs (681 up-regulated and 163 down-regulated) were screened out (Figures 4(d) and 4(e)). e detail information of each DEG was detailed in Table S2. We subsequently conducted CMap via these DEGs and further screened out seven molecule drugs (Table 2). Among them, sulpiride (P � 0.00008), famotidine (P � 0.00008), and nalidixic acid (P � 0.0118) showed strong potential to treat PAAD.

Establishing an Immune-Related Prognostic Signature.
Univariate Cox analysis of OS of these four IRGs was preliminary conducted (Figure 10(a)). By using Coxph function in R package "survival", we conducted Schoenfeld individual test for investigating the proportional hazards assumption.  Figure 5: Validation of hub genes. Kaplan-Meier survival curve based on GEPIA database revealed that PAAD patients with higher expression of hub genes had a significantly shorter overall survival time (S100A6: A; S100A10: E; S100A16: I; SDC1: (m)) and disease free survival time (S100A6: B; S100A10: F; S100A16: J; SDC1: (n)). Expressions of S100A6 (c), S100A10 (g), S100A16 (k), and SDC1 (o) in PAAD were significantly higher than these in normal tissues based on TCGA-PAAD database (P < 0.05). High expression of S100A6 (d), S100A10 (h), S100A16 (l), and SDC1 (p) related to higher tumor stage.
e global Schoenfeld test showed no significance (P � 0.8454, Figure 10(b)). Also, each variable including S100A6 (P � 0.3039), S100A10 (P � 0.5016), S100A16 (P � 0.6219), and SDC1 (P � 0.6655) owned no significance (P > 0.05, Figure 10(b)). us, this Cox model was conformed to the proportional hazards assumption. After this, all the prognostic biomarkers including S100A6, S100A10, S100A16, and SDC1 were used to establish the risk signature. e risk score of a sample was evaluated using the following equation: Risk score � 0.248 × Exp S100A6 + 0.355 × Exp S100A10 + 0.404 × Exp S100A16 + 0.240 × Exp SDC1 . Table S3 contained the risk scores of all TCGA-PAAD samples. 177 PAADs were split into two groups (high-risk group (n � 88), low-risk group (n � 89)) by setting the median risk score as cut-off standards. Next process concluded that PAAD patients with low-risk scores owned longer OS (Figure 10(c), P � 0.0074). Furthermore, the ROC values of this risk system were shown in Figure 10 By visualizing the distribution of patients in the two different groups (Figures 10(e) and 10(f )), we found that patients in high-risk groups were more likely to die comparing with these in low-risk group. With the aim of validating the repeatability and applicability of this signature, we repeated what we did in TCGA data via MMSD data, Table S4 contained the RSs of all PAADs. We also classified PAADs into high-(n � 134) and low-risk group (n � 135). PAADs who had higher-risk scores owned obviously poorer OS (Figure 10  Not measured

A Nomogram with Clinical Utility Was Constructed
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 analysis.   . In total, 11-12 samples were analyzed for S100A6, S100A10, S100A16, and SDC1. Immunohistochemistry (IHC) staining was evaluated as high/medium/low staining or not detected. (b) Representative IHC stainings for S100A6. (c) Representative IHC stainings for S100A10. (d) Representative IHC stainings for S100A16. (e) Representative IHC stainings for SDC1. All hub genes showed a membrane and/or cytoplasmic staining in tumor cells.
univariable Cox analysis suggested, risk score (P < 0.001), age (P < 0.007), tumor grade (P < 0.014), and stage (P � 0.030) were significantly associated with OS ( Figure 10(k)). Subsequent multivariate Cox analysis confirmed that the risk score could predict the prognosis of PAAD patients as individual (Figure 10(l)). Schoenfeld individual test was further performed for investigating the proportional hazards assumption as previous did. e global Schoenfeld test showed no significance (P � 0.107, Figure 10(m)). Also, each variable including age (P � 0.6687), tumor grade (P � 0.2792), pathologic stage (P � 0.0706), and risk score (P � 0.3409) was not statistically significant (P > 0.05, Figure 10(m)). us, this Cox model was conformed to the proportional hazards assumption. With the aim of applying clinical practice of this signature, we further constructed a nomogram based on risk score and age (Figure 11(a)), which showed significance in multivariate Cox analysis (Figure 10(k)). As the calibrate curve concluded, our nomogram could effectively predicate the survival rate of PAAD patient (Figures 11(b)-11(d)), least of all from long-term mortality (3-year OS (Figure 11(c)); 5-year OS (Figure 11(d)). Clinical net benefit of this nomogram was measured by DCA. e results suggested that the risk-signature-based nomogram owned better net benefit comparing to the norisk-score nomogram when predicting 1-(0.10 < Pt < 0.30, Figure 11(e)), 3-(0.20 < Pt < 0.65, Figure 11(f )), 5-year (0.20 < Pt < 0.65, Figure 11(g)) survival probability. By drawing ROC curves, the nomogram was validated to show great potential for PAAD patients OS prediction (Cindex: 0.780; AUC: 0.909; Figure 11(h)). Subsequent analysis concluded that this nomogram with immunerelated prognostic index revealed great steadiness across five years.

Identification of the Risk Signature Associated KEGG
Signaling Pathways. We further conducted GSEA to screen the abilities of the established index. Based on the standards we set in methods part, the risk score played significant role in Homologous recombination (Figure 12(a)), Pentose and glucuronate interconversions (Figure 12(b)), Ascorbate and aldarate metabolism (Figure 12(c)), and Linoleic acid metabolism (Figure 12(d)). e detail information for these pathways were shown in Figures 12(a)-12(d).

Role of the Immune-Related Prognostic Signature in
Immunotherapeutic Benefit Prediction. Immunotherapies (including PD-L1, PD-1 blockade, and etc) has undoubtedly emerged a major breakthrough in tumor treatment. In IMvigor210 (Table S5 showed the RSs of patients from IMvigor210), PAAD patients classified into high-risk score group (n � 147) owned significantly shorter survival (P < 0.001, Figure 13(a)), compared with patients in low-risk score group (n � 147). e risk score was further proved to show good potential for anti-PD-L1 immunotherapy predication by using IMvigor210 (Figures 13(b)-13(d)). Patients in lower risk score group were more likely to benefit from anti-PD-L1 treatment (Figures 13(b) and 13(c)), which was validated by Kruskal-Wallis test (P � 0.0059, Figure 13(d)). Next ROC analysis (based on IMvigor210) was used to evaluate the efficacy of anti-PD-L1 treatment by risk score, which recommended that risk score was a predictive biomarker of immunotherapeutic benefits (AUC: 0.714, Figure 13(e)). In addition, we also tried to find some relationship between risk score and anti-PD-1 treatment via GSE78220 cohort. e results suggested that PAAD patients in low-risk group exhibited significantly clinical benefits and an obviously prolonged survival (P � 0.024, Figure 13(f)). Also, patients in lower risk score group were more likely to benefit from anti-PD-1 treatment (Figures 13(g) and 13(h)). ere was a trend

16
Journal of Oncology that patients with response (CR/PR) to anti-PD-1 treatment occupied lower risk score (P � 0.742, Figure 13(i)). But unfortunately, further ROC analysis also demonstrated that risk score might not be an appropriate predictive tool to anti-PD-1 treatment benefits (AUC: 0.500, Figure 13(j)).

Correlation of Immune-Related Prognostic Signature with Immune Infiltration Level, And Tumor Status in PAAD.
Hopkins et al. proved T cell receptor repertoire features might act as a novel biomarker for PAAD immunotherapy response [42]. us, we also tried to explore the association among our index and immune cells. According to Figure 13(k), the signature was correlated to B cell naïve, monocyte, CD4+ T cell memory activated and CD8+ T cell, negatively. Meanwhile positively correlated to macrophage M0, myeloid dendritic cell activated, and T cell regulatory (Tregs). e P values were shown in Figure 13(k). Moreover, because of the present situation that the uncertainty of prediction of Tcell receptor (TCR) might to immunotherapy in metastatic cancer patients [43], we also attempted to explore whether the immune-related signature could predict it was a metastatic tumor or not. As the results showed ( Figure 13(l)), metastatic PAAD occupied higher-risk score compared with primary PAAD. But the result must be validated by using datasets with larger samples because the P value was greater than 0.05. Further ROC analysis also demonstrated that the immune-related prognostic signature might not effectively distinguish metastatic PAAD from primary PAAD (AUC: 0.548, Figure 13(m)).

Discussion
Pancreatic adenocarcinoma carries the lowest survival rate of all major organ cancers, which is of dismal prognosis and high mortality rate [44]. Most PAAD patients have no symptoms even in the advanced stage. Surgical management is the only treatment possible to cure PAAD, but its results is not satisfactory by now [44]. e 5-year survival rate of PAAD patients after receiving complete surgical resection is still less than 25% [45]. PAAD, with an insidious and atypical clinical symptoms, is a malignant neoplasm of the digestive tract that is difficult to diagnose and treat. PAAD in its early stages is not associated with a high rate of surgical mortality, whereas cure rates are low. Obviously, there was of great need of effective therapy target and prognostic biomarkers. erefore, screening novel prognostic biomarkers for PAAD and providing some novel choices for PAAD therapy are of urgent need.   21  20  19  18  17  16  15  14  13  12  11   20  30  40  50  60  70  80  90

Journal of Oncology
Cancer immunotherapy is an important tumor therapy potion for prevention and treatment of tumors and has attracted tremendous interests [46]. Nowadays, more and more researchers focus on screening out novel prognostic biomarkers related to immune microenvironment. However, similar studies regarding PAAD remains scarce. us, the aiming of the present study is to identify some biomarkers associated with prognosis of PAAD and explore some molecule drugs with therapeutic effect for PAAD.
Angelo et al. firstly investigated the correlation between PAAD patient prognosis and expression of 521 immune system genes [47]. Twenty immune system genes were carried out, which might influence PAAD prognosis. But there was a lack of validation for the 20 genes, and no further breakthrough had been made. Wu et al. had established a 3 IRG-based index, which showed definitive predictability for PAAD [48]. But not as we expected, the AUC for this model was less than 0.75, which might not be an effective tool for survival prediction. Chen et al. developed a two DEIRGbased signature, which could be used as an independent tool for the prognostic prediction of PAAD and to provide potential novel immunotherapy targets [49]. e similar question was that the AUC not as high as we expected (OS: 0.736) [49]. To learn from advantages and avoid disadvantages of these studies, we attempted to explored an immune-related prognostic signature via several authoritative methods and strict thresholds, based on multiple datasets and databases.
For the first time, we conducted WGCNA by using IRGs collected from ImmPort database. Four hub genes including S100A6, S100A10, S100A16, and SDC1 were screened out from the genes in the key module (yellow module). In addition, 22 overlapped DEGs between DEGs identified by MMD and EDGs identified by TCGA-PAAD were also selected. It was coincidence that the four hub genes identified by WGCNA were all the genes in the 22 DEGs. us, we regarded the four IRGs as potential prognostic biomarkers in PAAD and further verified them at different levels. All the four IRGs were overexpressed in PAAD tissues, compared with normal tissues, suggested by TCGA-PAAD data, MMD, GEPIA, and Oncomine database. Previous studies also indicated that the four IRGs were higher expressed in some cancer types. S100A6 has been proved to be highly expressed in epithelial cells, fibroblasts, and some cancer cell types, function of which was still unclear [50,51]. Serum S100A6 level was significantly increased in some solid carcinomas, including gastric carcinoma [52], bladder carcinoma [53] and ovarian carcinoma [54]. S100A10 was highly expressed in gastric cancer compared with normal gastric mucosa tissues,   which was involved in the occurrence and development of gastric cancer [55]. S100A16 was overexpressed in lung cancer [56], and colorectal cancer [57], which might play role in promoting the proliferation and migration of tumors. SDC1 was a kind of heparan sulphate proteoglycan, which was an important cell surface adhesion molecule to maintain cell morphology [58]. As reported, the disorder of SDC1 expression could effectively influence tumor cells invasion and metastasis [59]. S100A6, S100A10, and S100A16 are members of the S100 family of proteins containing 2 EF-hand calcium-binding motifs [60]. Some protein types encoded by S100 family have been used as tumor markers in clinical [60]. In conclusion, this study suggested that the expressions of S100A6, S100A10, S100A16, and SDC1 were upregulated in PAAD, which might be crucial biomarkers in the progression of PAAD. Furthermore, we validated the prognostic value of the four IRGs by using TCGA-PAAD data and MMSD. e results conducted that higher expression of the four IRGs were related to worse survival (OS and DFS) of PAADs. To learn from advantages and avoid disadvantages of these studies, we attempted to explored an immune-related prognostic signature via several authoritative methods and strict thresholds, based on multiple datasets and databases. Previous studies also indicated that over expression of serum S100A6 level was closely associated with the occurrence, development, prognosis and treatment of tumors [52]. It also had been reported that the up-regulation of S100A10 in renal cell carcinoma and bladder carcinoma was closely related to the poor prognosis, which might be correlated with the promotion of tumor cell proliferation, migration and invasion by S100A10 [61]. e disorder of SDC1 expression was also significantly correlated to poor prognosis of cancers [62]. As a conclusion, the present study supported the view that the four IRGs were closely associated with prognosis of tumors and might be novel immune-related prognosis biomarkers in PAAD.
Moreover, we established an immune-related prognostic signature base on the four prognostic biomarkers. It was worth mention that our immune-related signature might be the first one constructing (by combing DEG and WGCNA) for prognosis of PAAD patient prediction. is signature was validated to perform as an independent prognostic index with excellent potential of prognosis prediction of PAAD patients. With the aim of making the risk signature to become a clinical reality, we constructed a nomogram relying on the prognostic signature and age. e risk-signature-based nomogram was immediately validated to serve as a predictor for OS possibility of PAADs.
Argentiero et al. revealed that immune treatment might be effectively strategy through WNT pathway [63]. ese findings hold the promise to identify novel immune-based therapeutic strategies targeting WNT to enhance PDAC cytotoxicity and restore anti-PDAC immunity in nodepositive disease. Russano et al. discussed the current and possibly future applications of blood-based liquid biopsy in oncology, its advantages and its limitations in clinical practice [43]. Moreover, they specifically revealed its role as a tool to capture tumor heterogeneity in metastatic cancer patients. In another study, Hopkins et al. proved T cell receptor repertoire features associated with survival in immunotherapy-treated PAAD [42]. Whether these repertoires might predict response to immunotherapy in metastatic cancer patients is still unclear. us, we tried to explore the role of our signature in immunotherapeutic benefit prediction.
e results obviously demonstrated that this signature was related to response to anti-PD-1/L1 immunotherapy.
Some small molecule drugs were also identified, which could be useful in the treatment of PAAD. 844 DEGs were firstly identified by using MMD developed by four  Journal of Oncology independent datasets. We then performed CMap analysis and the results demonstrated that three molecule drugs including sulpiride, famotidine, and nalidixic acid might be potential choices in the treatment of PAAD. Some limitations might exist in the present study. As a multiple-dataset-based study, though we tried our best to design this bioinformatic study reasonably and logically, this research was lack of external experimental verification. Some experiments must be conducted in our further research to measure the relayed mechanisms of hub genes for PAAD. Besides, we will also validate the therapy potential of the three drugs in the subsequent analysis. In addition, clinical trials based on the extensive genomic technologies made the findings more reliable.

Conclusions
All in all, hub genes in WGCNA was firstly screened. en DEGs were further identified by using TCGA-PAAD data and MMD. Furthermore, four IRGs were screened out by overlapping core genes identified by WGCNA and abnormally expressed IRGs. All the four IRGs were screened out and verified to be correlated to dismal prognosis of PAAD. Moreover, 3 molecule drugs (sulpiride, famotidine, and nalidixic acid) were chosen and validated to have potential to treat PAAD. In conclusion, the present study indicated that four IRGs might be novel immune-related prognostic biomarkers of PAAD. Besides, three molecule drugs might have worked in the treatment of PAAD. An immune-related prognostic signature was constructed and validated as an independent prognostic biomarker with excellent potential of PAAD prognosis prediction. e nomogram based on the risk signature could act as a visual tool for OS probability prediction of PAADs.
Data Availability e data that support the findings of this study are openly available in Gene Expression Omnibus (GEO) database at http://www.ncbi.nlm.nih.gov/geo/, and e Cancer Genome Atlas (TCGA) database at https://genomecancer.ucsc.edu/.

Conflicts of Interest
ere are no conflicts of interest in the present study.

Authors' Contributions
Y. Y. put forward the study thought, the main content, and research methods. X. Z. and X. L. completed the whole analytic process of this study. Y. Y., X. L., J. X., Q. Z., and X. Z. organized and presented the results. X. Z. wrote the original draft of this manuscript. All authors agreed to publish this article. Xiao Zhang and Xiaomian Li contributed equally to this work.