RNA-sequencing (RNA-seq) data and clinical follow-up information data from TCGA-PAAD samples were downloaded from the TCGA database. The expression data and clinical information from the GSE57495, GSE62452 and GSE28735 datasets were downloaded from the GEO database. A total of 97 invasion-related genes were collected from the CancerSEA website ( Supplement Table 1 ).

The RNA-seq data from the TCGA-PAAD dataset was processed through the following steps: 1) Samples with no clinical follow-up information were removed; 2) The ENSEMBL gene IDs were converted to the Gene Symbol format; 3) The median value was calculated with multiple Gene Symbol expressions.

The following steps were used to process the GEO dataset: 1) Samples without clinical follow-up information were removed; 2) The probe IDs were converted to the Gene Symbol format; 3) Probes that corresponded to multiple genes were removed. 4) When multiple probes correspond to one gene, take the average value as the gene expression.

After preprocessing, we enrolled 176 samples from TCGA-PAAD, 63 samples from GSE57495 data set, 66 samples from GSE62452 data set, and 42 samples from GSE28735 dataset. The clinical characteristics of the patient samples are listed in Table 1 .

The expression profiles of 97 invasion-related genes were extracted from the TCGA-PAAD dataset, and univariate Cox regression analysis was performed to select significant prognostic genes using coxph function in R (p < 0.05). Next, the genes with significant results from the univariate Cox analysis were clustered using ConsensusClusterPlus (V1.48.0; parameters: reps = 100, pitem = 0.8, pfeature = 1, and distance = “Canberra”). The Pam and Canberra distances were used as a clustering algorithm and distance measure, respectively.

The gene set c2.cp.kegg.v7.0.symbols.gmt was selected, and significantly enriched pathways between different molecular subtypes were analyzed by GSEA. PAAD samples were divided into either a C1 or C2 subtype based on gene expression data from the TCGA-PAAD dataset in the GSEA input file. The thresholds for pathway enrichment analysis were p < 0.05 and false discovery rate (FDR) < 0.25.

DEGs between C1 and C2 subtypes were calculated using the limma package (6), and the filtering thresholds were FDR < 0.05 and | log 2 fold-change (FC) | > 1. The identified DEGs were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) enrichment analysis using the WebGestaltR (v0.4.2) package in R software.

Univariate Cox regression analysis was performed for each DEG (538 in total) using the coxph function in R to identify prognostic genes, and p < 0.05 was selected as the threshold for filtering. Lasso regression analysis was performed to further reduce the number of genes in the risk model using the glmnet package in R (7). In stepwise regression analysis, the selection of the model starts with the most complex model from which one variable is removed at a time to reduce the number of parameters according to the Akaike Information Criterion (AIC). The smaller the p-value of the regression model, the more superior the model. This indicates that the regression model fits the data well with fewer parameters. The prognostic model is made fit for clinical applications by performing stepwise regression to further reduce the number of genes.

The prognostic model was constructed based on the following equation:

where n refers to the number of genes identified for the multivariate Cox regression model; exp(Gi) is the expression value of gene i; and βi is the coefficient for gene i.

To verify the expression of the candidate three genes, tissue microarrays (TMA) comprised of 120 PAAD tissues and 30 normal samples were obtained from Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). The clinicopathological details of 120 PAAD tissues were shown in Table 2 . The studies were conducted in accordance with the International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS), and the research protocols were approved by the Ethics Committee of Hainan General Hospital, Hainan Affiliated Hospital of Hainan Medical University.

The TMA slides were dried overnight at 37°C, dewaxed in xylene, and dehydrated in a gradient ethanol series. Antigens retrieval was performed by heating the tissue sections in a microwave oven inside a vessel filled with EDTA antigen retrieval buffer (pH 9.0). Subsequently, the tissue sections were immersed in 3% hydrogen peroxide for 25 min to block the activity of endogenous peroxides. Next, the TMA tissues were coated with 3% bovine serum albumin (BSA) and sealed at room temperature for 30 min to reduce non-specific staining. Then, the TMA slides were incubated with anti-LY6D (1: 200 dilution; Novus Biologicals, NBP1-84029), anti-BCAT1 (1:50 dilution; Abcam, ab197941), and anti-ITGB6 (1:10 dilution; Abcam, ab197672) overnight at 4°C.

The tissues were rinsed with 0.01 mol/L phosphate buffer saline (PBS; pH = 7.4) for 5 min each. The tissues were incubated at room temperature for 50 min with horseradish peroxidase (HRP)-labeled goat anti-rabbit secondary antibody (1:200 dilution, ServiceBio, GB23303). Then, the tissues were washed in PBS and stained with 3,3-diaminobenzidine (DAB). Finally, the TMA sections were counterstained with Mayer’s hematoxylin, dehydrated, and fixed. To evaluate IHC staining, semi-quantitative scoring criteria were used.

The stained sections were scored by three pathologists who were blinded to the patients’ clinical characteristics. The scoring system was based on the proportion of positively stained cells in all tissues and the staining intensity of these positively stained cells. The staining intensity was classified as follows: 0 (negative), 1 (weak), 2 (moderate), or 3 (strong). The staining ratio of positive cells was classified as follows: 0 (<5%), 1 (5%-25%), 2 (26%-50%), 3 (51%-75%), or 4 (> 75%). According to the staining intensity and the proportion of positively stained cells, the tissues were graded as follows: 0-1 grade, negative (-); > 1-4, weakly positive (+); > 4-8, moderately positive (++), and > 8- 12, strongly positive (+++).

Thirty-five genes were found to be significantly associated with the prognosis of pancreatic cancer using univariate Cox analysis ( Supplement Table 2 ). Consistent cluster analysis showed that the samples could be clustered together at k=2 ( Figures 1A, B ). The expression levels of the invasion-related genes were significantly different between the C1 and C2 subtypes, and most genes were overexpressed in the C1 subtype ( Figure 1C ). The relationship between the subtypes and prognosis was further analyzed, and results showed that there were significant differences in survival times between the C1 and C2 subtypes ( Figures 1D, E , log-rank p < 0.05).

The results of the GSEA analysis showed the activation of more tumor-related pathways in the C1 subtype, such as pathways in cancer, notch signaling pathway, focal adhesion, extracellular matrix (ECM)-receptor interaction, and TGF-β signaling pathway ( Figure 1F ), suggesting that the C1 subtype is more closely related to cancer than the C2 subtype.

According to the thresholds mentioned in the methods section, 538 DEGs were obtained, of which 531 genes were upregulated and 7 genes were downregulated ( Supplement Table 3 ). The results demonstrated that the C1 subtype contains more upregulated genes than the C2 subtype. The volcano map of upregulated and downregulated DEGs between the two subtypes is shown in Supplementary Figure 1A . The expression patterns of the top 50 upregulated DEGs and all the downregulated DEGs were shown in a heatmap ( Supplementary Figure 1B ). The results of the GO enrichment analysis of DEGs showed that 548 Biological Process (BP) terms were significantly different between the two subtypes (FDR < 0.05). The first 15 BP terms were plotted (FDR < 0.05), as shown in Supplementary Figure 1C . The first 15 Cellular Component (CC) terms were plotted, as shown in Supplementary Figure 1D . Fifty-two Molecular Function (MF) terms were significantly different between the two subtypes (FDR < 0.05). The results of the first 15 MF terms are shown in Supplementary Figure 1E . The KEGG pathway analysis of DEGs showed 27 significantly enriched pathways (FDR < 0.05). Further visualization of the top 10 enriched pathways showed that genes were significantly enriched in tumor-related pathways such as the ECM-receptor interaction pathway, focal adhesion, and the PI3K-Akt signaling pathway ( Supplementary Figure 1F ).

To identify the relationship between molecular subtypes and immune scores in the TCGA-PADD dataset, the ESTIMATE package was used to evaluate the three immune scores: stromal, immune, and estimate scores. MCPcounter was used to evaluate 10 types of immune cells, and the single-sample GSEA (ssGSEA) method in the GSCA package was used to evaluate 28 types of immune cells (8). Meanwhile, the difference in immune scores between the two molecular subtypes was compared. The results showed that the immune scores of the C1 subtype were higher than those of the C2 subtype ( Figures 2A – C ). The heatmap of the immune scores of the two subtypes is shown in Figure 2D .

By performing univariate Cox analysis of the DEGs between the C1 and C2 subtypes, 18 prognostic genes were identified. Lasso regression analysis was performed to further reduce the number of prognostic genes. The locus of each independent variable is shown in Supplementary Figure 2A . As the value of lambda (λ) increased, the number of independent variables tending to zero also increased. A 10-fold cross-validation was performed to construct the model, and the confidence interval under each λ is shown in Supplementary Figure 2B . The model was found to be optimal when λ = 0.05667557, so a λ of 0.0567557 was chosen for further analysis of the prognostic genes. Six genes, namely LY6D, DKK1, BICC1, BCAT1, ITGB6, and PTGES were identified as the hub genes when λ = 0.0567557. The number of model genes was further reduced by stepwise regression, and finally, three genes were obtained: LY6D, BCAT1, and ITGB6. The risk score based on the final 3-gene prognostic model was calculated as follows: Risk score = 0.1627483 * LY6D + 0.2210480 * BCAT1 + 0.2005339 * ITGB6.

Risk scores of each sample were calculated based on the expression level of LY6D, BCAT1, and ITGB6, and a risk score distribution was plotted for each sample, as shown in Figure 3A . The results showed that a higher risk score was associated with worse outcomes, and high expression levels of LY6D, BCAT1, and ITGB6 were associated with a higher risk score. The timeROC package was used to analyze the receiver operating characteristic (ROC) curve of risk score; the 1-, 2-, and 3-year predictive classification efficiencies were 0.76, 0.78, and 0.75, respectively, as shown in Figure 3B . The samples were divided into a high-risk group and a low-risk group based on the risk scores. Finally, 50 and 56 samples were placed into the high- and low-risk groups, respectively. The KM curve showed a significant difference in the expression of DEGs between the high- and low-risk groups (p < 0.01) ( Figure 3C ).

Further analysis of the relationship between the risk score and clinical features showed that the 3-gene prognostic model could significantly distinguish between age, sex, TNM stage, clinical stage, tumor grade, alcohol consumption, chemotherapy, and radiation therapy between the high- and low-risk groups ( Figures 7A – P , p < 0.05). This suggests that the model also has good predictive power in distinguishing different clinical features.

The distribution of the risk score among the clinical features of the two groups was further compared. The results demonstrated that the risk score is significantly different between the N-stage and tumor grade (p < 0.05). The higher the tumor grade, the higher the risk score ( Figure 7Q ). The risk score of N1 was significantly higher than that of N0 ( Figure 7R ). The risk score of the C1 subtype with a poor prognosis was significantly higher than that of the C2 subtype with a good prognosis ( Figure 7T ). Moreover, the risk score was significantly different among existing immune molecular subtypes ( Figure 7S ).

In the TCGA-PAAD dataset, the univariate Cox regression analysis showed a significant correlation between the risk type and survival, while the multivariate Cox regression analysis showed a significant correlation between the risk score (Hazard ratio [HR] = 1.94, 95% confidence interval [CI] = 1.25–3.01, and p = 0.003) and survival. These results demonstrate the good predictive performance of the identified 3-gene prognostic model in clinical applications. Furthermore, the N stage (HR = 4.17, 95% 1.16–14.93, and p = 0.028) and grade (HR = 3.06, 95% CI = 1.3–7.21, and p = 0.011) were identified as independent prognostic risk factors for patients with pancreatic cancer. Chemotherapy (HR = 0.13, 5% CI = 0.05–0.36, and p < 0.001) was identified as an independent prognostic protective factor ( Figures 8A, B ).

The nomogram, which displays the results of the risk model directly and effectively, can be conveniently applied to the prediction of an outcome. The nomogram uses the length of the line to indicate the degree of influence that different variables have on the result and the influence of different values of variables on the result. According to the results of the univariate and multivariate analyses, the nomogram was constructed with the following clinical features: N stage, tumor grade, chemotherapy, and risk score ( Figure 8C ). The results showed that the risk score has the greatest effect on survival prediction, indicating that the risk model based on the 3-gene signature can accurately predict the prognosis of patients with pancreatic cancer. A calibration diagram was used to visualize the nomogram. The results showed that the nomogram performed well in determining the prognostic risk of patients with pancreatic cancer ( Figure 8D ).

Four prognostic risk models, including 15-gene signature (Chen) (9), 7-gene signature (Cheng) (10), and 6-gene signature (Stratford) (11) models, were compared with the identified 3-gene prognostic model. To facilitate comparison among the models, the risk score of each TCGA-PAAD sample was calculated using the same method, and the risk score was zscored according to the corresponding gene in all three models. Genes with a risk score greater than zero were categorized into a high-risk group and those with a risk score less than zero were categorized into a low-risk group. The prognosis difference between the two groups was further analyzed. There were significant differences in outcomes between the high-risk and low-risk groups in all three risk models ( Figures 9B, D, F , log-rank p < 0.05), the area under the curve (AUC)s at 1-, 2-, and 3-year of Cheng and Stratford models were lower than that of our model ( Figures 9C, E ). Although our 1-year AUC is smaller than the Chen model (0.73 vs 0.74), the AUC at 2 and 3 years is larger than his (0.77 vs 0.76, 0.81 vs 0.78, respectively) ( Figure 9A ) ( Supplement Table 4 ). Therefore, the 3-gene signature identified in this study represents a more reasonable and efficient model to determine the prognostic risk of patients with pancreatic cancer with the use of fewer genes.

The differences in the expression of the LY6D, BCAT1, and ITGB6 genes in PAAD and adjacent tissues were investigated. The expressions of LY6D, BCAT1, and ITGB6 in 120 cases of pancreatic cancer and 30 cases of para-carcinoma were detected by immunohistochemistry. The results showed that BCAT1, LY6D, and ITGB6 were significantly overexpressed in cancer tissues ( Figures 10A – C ). Many cases in the TMA cohort were not effectively followed up. Therefore, to compensate for this limitation, the Kaplan-Meier plotter database was used to obtain 177 samples with overall survival data and 69 cases with recurrence-free survival data. The results showed that patients with high expression of LY6D, BCAT1, and ITGB6 genes have a significantly worse prognosis than those with a low expression both in terms of overall survival and recurrence-free survival ( Figures 10D – I ). Our immunohistochemical results demonstrated that LY6D, BCAT1, and ITGB6 proteins were all overexpressed in PAAD samples compared to normal samples. Therefore, it can be speculated that these genes act as oncogenes in pancreatic cancer, and the upregulation of these genes is associated with a significantly worse prognosis in patients with pancreatic cancer.

A flowchart has been drawn to allow readers to better understand the research process of this study ( Figure 11 ).