Patient inclusion criteria in the TCGA dataset were [1] PAAD patients with a surgical operation of R0 resection [2], patients with no new tumor event or with liver metastasis during follow-up. The exclusion criteria were [1] patients with survival period <30 days or with incomplete survival data [2], patients with missing data of gene expression. Sixty-five patients were included in this study. During the follow-up period, the numbers of patients with liver metastasis and with no new tumor event were 15 and 50, respectively. The corresponding clinicopathological information of the PAAD patients was obtained via Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). The mRNA expression data [Fragments Per Kilobase of exon model per Million mapped fragments (FPKM)] of the PAAD samples in the TCGA dataset and the normal pancreas samples in the GTEx dataset were reprocessed using UCSC xena (https://xenabrowser.net/) to avoid data imbalance and transformed into log₂(FPKM+1). The batch effects between the TCGA and the GTEx cohorts were minimized by “limma” R package. To assess the batch effects, the variation of housekeeping gene expression between the TCGA and the GTEx datasets was evaluated as Eisenberg et al. reported [15]. The mRNA expression profiles of the 65 PAAD and 165 normal pancreas samples were obtained for downstream analysis. The GRCh38 file was downloaded from the Ensembl website (https://asia.ensembl.org/index.html) for annotation of the mRNAs.

The 65 PAAD patients in the TCGA cohort were training dataset. The whole TCGA-PAAD dataset with 147 patients (cases with a survival period ≥30 days and complete gene expression data), GEO (https://www.ncbi.nlm.nih.gov/geo/) and ICGC (obtained from SurvExpress website, http://bioinformatica.mty.itesm.mx:8080/Biomatec/SurvivaX.jsp) were used as validation datasets.

Based on the threshold of |logFC| > 1 and p value <0.05, the differentially expressed genes (DEGs) of normal pancreas samples vs. primary tumor samples, and DEGs of primary tumor samples without new tumor event vs. primary tumor samples with liver metastasis during follow-up were screened by Wilcoxon test. Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify the intersection of DEGs in PAAD tissues and PAAD tissues with liver metastasis.

Prognostic Gene Signature Construction Based on Liver-metastasis-related mRNAs.

Univariate regression Cox analysis was performed to filter mRNAs associated with OS. Then the identified mRNAs were included in the least absolute shrinkage and selection operator (LASSO) Cox regression model to develop a liver-metastasis-related prognostic signature (LM-PS) for the PAAD patients involving seven liver-metastasis-related genes and derive the regression coefficient of each gene by “glmnet” and “survival” R packages. Thereafter, risk scores for patients were calculated according to the following formula: R i s k   s c o r e = ∑ i = 1 n ( C o e f i ∗ E x p i ) where E x p i represents the expression of each prognostic mRNA in the LM-PS, and C o e f i represents the coefficient of the corresponding mRNA. PAAD patients were grouped into low- and high-risk groups based on the median value of the risk scores. Principal component analysis (PCA) was performed for the two groups using the R package “factoextra”.

In the TCGA dataset, the liver-metastasis-related prognostic signature and corresponding clinicopathological factors of the PAAD patients were included in univariate/multivariate Cox proportional hazards analysis and logistic regression model analysis to identify the independent predictors of OS and liver metastasis, respectively. A nomogram was established based on the independent predictors of OS with the “rms” R package. The receiver operating characteristic (ROC) curves and the area under the curve (AUC) values were used to evaluate the predictive ability of the established models in this study.

Based on the threshold of |logFC| > 1 and p value <0.05, 408 DEGs between high- and low-risk groups were identified and uploaded to the “Metascape” website (https://metascape.org/) for enrichment analysis involving Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG).

The infiltration levels of the 22 immune cells were determined to assess the tumor microenvironment between high- and low-risk groups and between liver-metastasis and non-new-tumor groups in the TGCA dataset by the CIBERSORT R script v1.03. A p value <0.05 was set as the significance threshold.

The cBioPortal website (http://www.cbioportal.org/) is a large depository that provides access to cancer genomics data carried out by many institutions including the TCGA database [16]. In this study, the cBioPortal was used to explore the connection between the gene mutation and the prognosis of the 65 PAAD patients in the TCGA dataset.

The seven liver-metastasis-related genes identified by the Lasso Cox regression model were uploaded to NetworkAnalyst (https://www.networkanalyst.ca/faces/home.xhtml) to construct a transcription factor‐microRNA coregulated network. The transcription factor‐microRNA coregulatory data were predicted from the Regulatory Network Repository. The transcription factor‐microRNA coregulatory network was visualized by Cytoscape 3.8.1.

PAAD patients in the TCGA cohort were divided into high- and low-risk groups based on the median value of the risk scores. Mean ± standard deviation or medians (with interquartile range), and counts with percentages were used to present continuous variables and categorical variables, respectively. Continuous variables were compared by Student’s t-test or Mann-Whitney U-test, and categorical variables by Chi-square test. Kaplan-Meier curves were plotted to analyze patient survival. The log-rank tests were performed to analyze the differences in patient survival. The statistical analysis was carried out with the R software 4.0.2, SPSS 22.0, and GraphPad Prism 8.0.1. p values <0.05 were considered to indicate statistical significance (two-sided).

This study included 65 PAAD patients with R0 resection from the TCGA dataset, and the corresponding clinicopathological factors of the patients were presented in Table 1. As shown in the table, with the exception of survival status and pathologic M stage, there was no significant difference in terms of all the other clinicopathological variables that were analyzed. Previous studies successfully compared the gene expression between the GTEx and TCGA datasets, and minimized the batch effects in data processing [14,17]. To assess the batch effects in this study, the variation of housekeeping genes expression between the TCGA and GTEx datasets was compared, and a high correlation between the PAAD and normal pancreas tissues was observed (Pearson R = 0.89, p < 0.0001, Supplementary Figure S1), indicating minimal batch effects between the two datasets. Subsequently, we screened DEGs based on the standards of |logFC| > 1 and p < 0.05 between the 65 PAAD tissues vs. 165 normal pancreas tissues and between the 15 primary PAAD tissues with liver metastasis vs. 50 primary PAAD tissues without new tumor event. The results showed that 4,996 and 62 genes were differently expressed in PAAD tissues and PAAD tissues with liver metastasis, respectively. Thereafter, the Venn diagram web tool was used to get the intersection of the two groups of DEGs, and 45 liver-metastasis-related genes were identified. The study flow chart is shown in Figure 1.

To build the LM-PS for forecasting the OS of PAAD patients, we included the expression of the 45 liver-metastasis-related genes into the univariate Cox regression analysis and identified 33 genes significantly associated with OS (Supplementary Table S1). Subsequently, based on the 33 prognostic genes, the LASSO Cox analysis (Figures 2A,B) was performed to construct a seven-gene LM-PS containing anoctamin-1 (ANO1), family with sequence similarity 83, member A (FAM83A), G protein-coupled receptor 87 (GPR87), integrin beta-6 (ITGB6), kallikrein-10 (KLK10), serine protease inhibitor, clade E member 1 (SERPINE1), and small integral membrane protein 32 (SMIM32). The detailed information and the prognostic ability of the seven liver-metastasis-related genes were presented in Table 2 and Figure 2C, respectively. Next, the risk score for predicting OS was calculated for each patient in the training dataset based on the expression and coefficients (Figure 2D) of the seven liver-metastasis-related genes. Then the patients were separated into low- and high-risk subgroups based on the median value of the risk scores (Figure 2E). The survival status of the patients in the high- and low-risk groups was visualized in Figure 2F showing that the mortality rate post R0 resection was higher in the high-risk group than in the low-risk group. The heatmap showed that the expression of the risk genes (ANO1, FAM83A, GPR87, ITGB6, KLK10, and SERPINE1) were up-regulated, while the expression of the protective gene (SMIM32) was down-regulated (Figure 2G) with the increasing risk score. In addition, the PCA was used to assess the biological difference between low- and high-risk groups based on the expression of the seven liver-metastasis-related genes. The result demonstrated that the low- and high-risk patients were distributed separately in distinct directions (Figure 2H).

The patients in the training dataset were divided into two groups based on the median expression of each liver-metastasis-related gene. The Kaplan-Meier curves showed that high expression of ANO1, FAM83A, GPR87, ITGB6, KLK10, SERPINE1, and low expression of SMIM32 were significantly associated with worse patient survival (Figures 3A–G). In addition, based on the two different risk groups distinguished by the LM-PS, the Kaplan-Meier curve also showed that patients with high-risk had a significantly lower OS rate than those with low-risk (p = 6.093E-05, Figure 3H). The OS predictive ability of the LM-PS was confirmed by the ROC curves in the TCGA cohort with 1-year AUC = 0.864, 2-years AUC = 0.899, and 3-years AUC = 0.856 (Figure 3I).

The distribution of clinicopathological factors with the risk score increasing is shown in Figure 3J. In addition, risk scores were observed to be higher in patients with G2 and G3/G4, pathologic stage II, and liver metastasis compared with the corresponding control groups but were not significantly associated with age, gender, pathologic T stage, or pathologic N stage (Figures 3K–Q).

As expected, the Kaplan-Meier curve showed that the OS rate was remarkably lower in PAAD patients with liver metastasis after R0 resection than those without new tumor events (Figure 4A). This study tried to provide an insight into the factors associated with the liver metastasis of PAAD. The risk score, age, gender, histologic grade, pathologic stage, pathologic T stage, pathologic N stage, and the number of positive lymph nodes were included in univariate logistic regression model analysis. Interestingly, only risk score was significantly associated with liver metastasis [Odds Ratio = 1.176, 95% confidence interval (CI) = 1.061–1.302, p = 0.002; Table 3]. The ROC curve confirmed that the risk score harbored a good ability to predict liver metastasis of PAAD in the TCGA cohort (AUC = 0.756, Figure 4B).

In the whole TCGA-PAAD cohort, the forest plot showed that high expression of ANO1, FAM83A, GPR87, ITGB6, KLK10, SERPINE1, and low expression of SMIM32 were significantly associated with worse patient survival (Figure 5A). Based on the LM-PS, risk scores were calculated for patients in the whole TCGA-PAAD and GSE57495 cohort and obtained from the SurvExpress website for the ICGC cohort, respectively. The clinicopathological data of the GSE57495 were presented in Supplementary Table S2. Patients in the three datasets were divided into high- and low-risk groups based on the median value of risk scores, respectively. The Kaplan-Meier curves showed that the high-risk patients had significantly lower OS rates than the low-risk patients in the whole TCGA-PAAD, GSE57495, and ICGC datasets (Figures 5B–D). Risk scores and survival status distributions of GSE57495 and ICGC were presented in Supplementary Figures S2A–D. Then the ROC curves were plotted to assess the prognostic values of the risk scores. The 1-, 2-, and 3-years AUCs were 0.766, 0.786, and 0.748 in the whole TCGA dataset, respectively (Figure 5E); 0.601, 0.710, and 0.645 in the GSE57495 dataset, respectively (Supplementary Figure S2E); and 0.719, 0.664, and 0.601 in the ICGC dataset, respectively (Supplementary Figure S2F). These results indicated that the LM-PS harbored a moderate OS predictive ability for PAAD patients.

Univariate and multivariate Cox regression analyses were used to assess the prognostic value of the LM-PS in the TCGA-PAAD dataset. In the training set, the risk score and clinicopathological factors of PAAD patients were included in univariate Cox analysis, and the results showed that risk score was significantly associated with OS [hazard ratio (HR) = 1.132, 95% CI = 1.077–1.190, p < 0.001; Table 4]. Furthermore, the multivariate Cox analysis demonstrated that risk score was an independent predictor of OS (HR = 1.126, 95% CI: 1.060–1.197, p < 0.001; Table 4).

Univariate and multivariate Cox regression analyses were also performed in the whole TCGA-PAAD dataset. The univariate Cox regression analysis showed that risk score was significantly associated with OS (HR = 2.668, 95% CI = 1.761–4.042, p < 0.001; Table 5). The multivariate Cox regression analysis further indicated that risk score was an independent predictor of OS in the whole TCGA-PAAD dataset (HR = 2.762, 95% CI: 1.680–4.541, p < 0.001; Table 5).

In the training cohort, the risk score, number of positive lymph nodes, and histologic grade were included to establish a clinically applicable quantitative nomogram to predict the OS of PAAD patients (Figure 6A). The calibration plot revealed that the predicted 1-, 2- and 3-years OS using nomogram agreed with the observed OS in the TCGA dataset (Figure 6B). Furthermore, ROC curves were used to assess the prediction accuracy of the nomogram both in the training cohort and the whole TCGA-PAAD cohort. The AUCs values for 1-, 2- and 3-years in PAAD were 0.947, 0.913, and 0.940 in the training cohort, respectively (Figure 6C); and 0.815, 0.826, and 0.791 in the whole TCGA-PAAD cohort, respectively (Figure 6D), indicating high accuracy of the nomogram in the TCGA-PAAD cohort.

Then the ROC curves were also plotted to assess the prognostic values of the nomogram and other independent predictors (risk score, histologic grade, and the number of positive lymph nodes) in the whole TCGA-PAAD dataset. The AUCs of the risk score, histologic grade, number of positive lymph nodes, and the nomogram were 0.753, 0.619, 0.665, and 0.815, respectively, indicating that the inclusion of risk score added to the accuracy of the other two factors included in the nomogram (Figure 6E).

The landscape of immune cell infiltration of PAAD in the TCGA cohort was explored between low- and high-risk groups and between liver-metastasis and no-new-tumor groups using the CIBERSORT algorithm. The violin plots showed that patients in the high-risk subgroup harbored significantly higher prevalence in M0 macrophages (p = 0.004) and Eosinophils (p = 0.003), but a lower prevalence in naive B cells (p = 0.022) than those in the low-risk subgroup (Supplementary Figure S3A). However, no immune cell proportions were observed to be different between non-new-tumor and liver-metastasis subgroups (Supplementary Figure S3B ).

The transcription factor‐microRNA coregulatory network based on ANO1, FAM83A, GPR87, ITGB6, KLK10, SERPINE1, and SMIM32 was obtained from NetworkAnalyst. This network contains the seven genes of the LM-PS, 63 transcription factors, and 37 microRNAs, with 86 associations between the LM-PS genes and transcription factors, and 38 associations between the LM-PS genes and microRNAs. As shown in Figure 7, the transcription factors FOS, JUND, USF1, and MYC were observed to regulate three LM-PS genes, the microRNA hsa-miR-224 was observed to modulate two LM-PS genes, and SERPINE1 was the most frequently modulated LM-PS gene. The increase of ANO1, FAM83A, GPR87, ITGB6, and the decrease of SMIM32 were significantly associated with the increase of SERPINE1, while no significant correlation was observed between the expression of KLK10 and SERPINE1 in the TCGA dataset (Supplementary Figures S4A–F).

We explored the genomic alteration for the TCGA-PAAD cohort in this study (data of 64 out of the 65 patients could be found on the cBioPortal website). In this study, 45.3% of the patients (29/64) had at least one LM-PS gene alteration with FAM83A, SERPINE1, and GPR87 being the most frequently altered genes (22%, 14%, and 11%, respectively; Supplementary Figure S5A). The oncoprints of low- and high-risk subgroups are shown in Figures 8A,B. Patients in the high-risk subgroup harbored significantly higher genetic alteration frequency in FAM83A (p = 0.002), GPR87 (p = 0.037), but lower frequency in SMIM32 (p = 0.024) than those in the low-risk subgroup (Supplementary Table S3). In total, no significant difference was observed in OS between the altered and non-altered subgroups (Supplementary Figure S5B). For the individual LM-PS gene alteration, compared to patients with no genetic alteration, patients with alteration of FAM83A and ITGB6 were significantly associated with lower OS (Supplementary Figures S5C,D), while those with alteration of SMIM32 were significantly associated with higher OS (Supplementary Figure S5E). The most frequent mutated genes between the patients with no LM-PS gene altered and those with at least one LM-PS gene altered were also explored. The results showed that the altered subgroup harbored higher enrichment of the V-Ki-ras2 Kirsten Rat Sarcoma Viral Oncogene Homolog (KRAS) (71.43% vs. 48.57%, p = 0.0571), Tumor Protein 53 (TP53) (64.29%% vs. 42.86%, p = 0.0748), and Cyclin-dependent Kinase Inhibitor 2A (CDKN2A) (28.57% vs. 8.57%, p = 0.0405) than the non-altered subgroup (Supplementary Figure S5F).

GO and KEGG enrichment analyses were performed for the 408 DEGs between high- and low-risk subgroups to identify the most relevant biological processes and pathways. The GO annotations showed that the DEGs enriched in GO terms mainly related to ECM, cell adhesion, and locomotion such as ECM organization, positive regulation of locomotion, cell-substrate adhesion, cell adhesion molecule binding, cell junction organization, collagen binding, regulation of cell adhesion, and ECM binding (Figure 9A). The most enriched KEGG pathway was also ECM-receptor interaction. Additionally, the KEGG enrichment analysis showed that DEGs were enriched in several tumor-related pathways such as proteoglycans in cancer, pathways in cancer, transcriptional misregulation in cancer, prostate cancer, and p53 signaling pathway (Figure 9B).