The transcriptome profiles of patients with PDAC were obtained from TCGA database (https://cancergenome.nih.gov/, accessed at April 20, 2017), and normalized using the DESeq Package in R (36,37). The corresponding clinical data were acquired from the University of California, Santa Cruz Xena (UCSC Xena; http://xena.ucsc.edu/, accessed at April 20, 2017). In order to eliminate interference from unrelated factors, the patients were selected based on the following inclusion criteria: i) Histological validation; ii) pathological stage I or II according to the 7th American Joint Committee on Cancer (AJCC); iii) availability of complete survival data; and iv) having undergone pancreaticoduodenectomy.

The CCR genes were functionally annotated using Gene Ontology (GO) terms, and the associated pathways were determined by Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/home.jsp version 6.8, accessed at January 19, 2018) (38). The correlation between CCR genes was analyzed by Pearson's correlation coefficient using the corrplot package in R (version 1.2.1335; www.r-project.org). A protein-protein interaction (PPI) map was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (https://string-db.org/, version 11, accessed at march 30, 2019) (39). All CCR gene symbols (CCR1-CCR10) were entered into the platform for Homo sapiens and an interaction score >0.4 was considered to be significant. Finally, a gene-gene interaction network was constructed using GeneMANIA (http://genemania.org/, accessed at March 30, 2019) (40).

The association between different clinical factors and the prognosis of patients with early-stage PDAC was determined using Kaplan-Meier analysis with the log-rank test; the relevant factors were then included in the multivariate Cox proportional risk regression model to identify the CCR genes significantly associated with OS. Based on the results of the survival analysis of individual genes, combined effect survival analysis was performed and a nomogram was constructed. The CCR genes associated with PDAC prognosis were assessed using combined effect survival analysis (Kaplan-Meier analysis with log-rank test) and the multivariate Cox proportional risk regression model. The nomogram was constructed in R (version 3.5.2; www.r-project.org) using the rms package, based on clinical variables and the expression levels of CCR genes. The scale marked on the line indicates the value range of each variable, and the length of the line segment reflects the contribution of this factor to the outcome event.

According to the results of the survival analysis, the CCR genes associated with PDAC prognosis were combined to construct a prognostic model based on gene expression level. The risk score formula was as follows: Risk score=expression of gene₁ × β₁ + expression of gene₂ × β₂+… expression of Gene_n × β_n (41,42), where β_n is the regression coefficient derived from the result of multivariate Cox proportional hazards regression analysis for the corresponding gene. Based on the median risk score value, the patients were divided into a high- and low risk group. To assess the predictive value of the prognostic signature, a time-dependent ROC curve was constructed using the survivalROC package in R (43). Survival analysis was performed to compare prognoses between the high- and low-risk groups.

In order to identify the pathways in which the CCR gene is enriched, and to determine whether the CCR genes in each gene set are enriched in the upper or lower part of the phenotype-related sorted gene list, genome-wide expression profile datasets and corresponding grouping files determined by the expression of CCR genes were uploaded to GSEA (44) for enrichment analysis with database c2 and c5 of the Molecular Signatures Database (MSigDB) (45). A set of genes with both false discovery rate (FDR) <0.25 and P<0.05 was considered to be statistically significant.

Statistical analysis was conducted using SPSS 22.0 (IBM Corporation) or R 3.52 (https://www.r-project.org/). Hazard ratios (HRs) and 95% confidence intervals (CI) were used to indicate the relative risk between the high-C and low-CCR expression groups. Pearson's correlation coefficient was used to determine the correlation between CCR genes, where P<0.05 was considered to indicate a statistically significant result. The FDR control in GSEA was achieved using the Benjamini-Hochberg procedure and adjusted for multiple testing (46–48), where an FDR<0.25 was considered to indicate a statistically significant difference.

The expression profiles of patients with PDAC were acquired from TCGA database. After screening based on the inclusion criteria, patients that fell outside of these parameters were eliminated, and the profiles of the remaining 112 patients were further analyzed.

As shown in Fig. 1A and Table SI, the results of GO and KEGG pathway analysis indicated that CCR genes are primarily involved in pathways related to immunity and inflammation, and that the JAK/STAT, mitogen-activated protein kinase (MAPK) and nuclear factor (NF)-κB signaling pathways were associated with chemokine signaling (Fig. S1). PPI and gene-gene interaction network analyses revealed that CCRs interact closely with each other (Fig. 1B and C). A matrix graph of Pearson's correlation analysis indicated that CCR1, CCR2, CCR4, CCR5 and CCR8 are closely related to each other with a correlation coefficient ≥0.5. In addition, a higher degree of correlation was observed among CCR4, CCR6 and CCR7 (correlation coefficient ≥0.7) (Fig. 2). The numbers in each grid represent the correlations between the corresponding genes.

The clinical data of all patients are summarized in Table SII. Histological grade, targeted molecular therapy, radiation therapy and residual resection were all significantly associated with OS. Patients harboring tumors of a higher histological grade, and those who did not receive either targeted or radiation therapies or undergo residual resection were at a higher risk of poor prognosis. In addition, high expression levels of CCR5 (adjusted P=0.012; adjusted HR=0.478, 95% CI=0.269–0.852), CCR6 (adjusted P=0.026; adjusted HR=0.527, 95% CI=0.299–0.927) and CCR9 (adjusted P=0.001; adjusted HR=0.374, 95% CI=0.209–0.670) were significantly associated with lower mortality rates (Table I and Fig. 3). The nomogram also indicated that CCR5, CCR6 and CCR9 may contribute to the prognosis of PDAC, with low expression corresponding to a high point (Fig. 4E).

The patients were stratified into groups based on the expression levels of different CCR genes. The expression levels of CCR genes in different groups are summarized in Table II. Favorable overall survival was observed in Group D (compared with Groups A, B and C; adjusted P=0.012; adjusted HR=0.434, 95% CI= 0.226–0.833), Group IV (compared with Group I, II and III; adjusted P<0.001; adjusted HR=0.236, 95% CI=0.107–0.520), Group d (compared with Group a, b and c; adjusted P=0.001; adjusted HR=0.284, 95% CI=0.136–0.595) and Group 4 (compared with Group 1, 2 and 3; adjusted P=0.001; adjusted HR=0.253, 95% CI=0.112–0.574) (Table II and Fig. 4A-D).

Using both single gene survival analysis and combined effect survival analysis, CCR5, CCR6 and CCR9 were demonstrated to be associated the prognosis of patients with early-stage PDAC; therefore, these three genes were selected for the construction of the prognostic signature. The regression coefficient of CCR5, CCR6 and CCR9 from the multivariate Cox proportional hazards regression model was −0.836, −0.618 and −0.476 respectively. Because all β-values in this investigation were <0, and to make the result easier to interpret, a constant was added to the end of the following risk score formula: Risk score=expression of CCR5 × −0.836 + expression of CCR × −0.618 + expression of CCR × −0.476 + 4[constant]. The effect of the constant is to ensure that the risk score output is >0. Survival analysis between the high and low risk score groups indicated that a high risk score was significantly associated with the poor outcome of patients with early-stage PDAC (adjusted P=0.018; adjusted HR=1.988, 95% CI=1.125–3.513) (Fig. 4F and G). Time-dependent ROC analysis demonstrated that the prognostic signature effectively predicted the outcome of patients with early-stage PDAC (1-year AUC=0.674; 2-year AUC=0.649; 3-year AUC=0.673; Fig. 4H).

Since CCR5, CCR6 and CCR9 were favorably associated with OS, the patients were stratified according to their respective median expression values. GSEA results are displayed in Figs. 5–7. Analysis of the C2 (curated) gene sets revealed that CCR5 was enriched in TP53 target, TP63 target, MAPK signaling pathway, generic transcription pathway, DNA damage and STAT5A target (Table SIII, Fig. 5A-F); in the C5 (GO) gene sets, CCR5 was enriched in inflammatory response, STAT cascade, MAPK cascade, regulation of NF-κB and endothelial proliferation (Table SIII, Fig. 5G-P); CCR6 was enriched in TGF-β1 signaling pathway, TP53 and TP63 targets, KEGG MAPK signaling pathway, MAPK14 targets and NF-κB signaling in the C2 set (Table SIV and Fig. 6A-G), and in STAT cascade, MAPK cascade, NF-κB import into nucleus, NF-κB signaling and transcription factor importing into nucleus in the C5 set (Table SV, Fig. 6H-M); CCR9 was enriched in the IL-2/STAT5 pathway, proliferation, NF-κB atypical pathway, STAT 5 targets and PTEN pathway in the C2 set (Table SVI and Fig. 7A-G).