A total of 12 PAAD tissues and the corresponding paracancerous tissues were collected from the Affiliated Hospital of Guizhou Medical University from August 2020 to April 2021. This study was approved by the Ethics Committee of Guizhou Medical University, and all participants with PAAD provided informed consent forms.

PAAD cell lines PANC-1, MIA PaCa-2, BxPC-3, AsPC-1, and SW1990 were purchased from the American Type Culture Collection (Manassas, VA, United States), and the normal human pancreatic duct epithelial (HPDE) cell line was purchased from Beijing North Carolina Chuanglian Biotechnology Research Institute (Beijing, China). All cell lines were cultured in the corresponding medium with 10% fetal bovine serum (Gibco, Waltham, MA, United States).

Gene expression profiles of 182 PAAD patients (178 tumor samples and 4 normal samples) were downloaded from TCGA (https://portal.gdc.cancer.gov/) database, including lncRNAs, mRNAs, miRNAs, HTseq (a python package) counts, fragments per kilobase of exon per million mapped fragments (FPKM), and the clinical information of each patient. Then, the PAAD group and the normal group were filtered through Perl language, the differentially expressed lncRNAs, mRNAs, and miRNAs were identified using the “DEseq2” package in R language (version 4.0.2), and heatmaps and volcano maps were plotted. DEG’s filtering criteria were as follows: fold change (FC) of differential expression of lncRNAs and mRNAs was logFC > 2, −logFC < −1 (miRNA) and p value < 0.01.

The correlation analysis of differentially expressed lncRNAs and mRNAs in TCGA database was performed with the “GDCRNATools” package in R language. Then, the differentially expressed miRNAs that can regulate lncRNAs and mRNAs were identified using the hypergeometric test and correlation analysis to construct a ceRNA network (p < 0.05 as the filter threshold criteria). The PAAD ceRNA network was visualized by introducing edge and nodal gene information into the interaction network with Cytoscape v3.8.0. Subsequently, differentially expressed genes (DEGs) in the ceRNA network were annotated with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The above analysis was performed with the “clusterProfiler,” “org.Hs.eg.db,” “enrichplot,” and “ggplot2” packages in R language, and then Cytoscape was employed for the GO/KEGG network analysis. Among them, GO analysis was utilized for identifying the biological characteristics, including biological processes (BPs), cellular components (CCs), and molecular functions (MFs). Meanwhile, KEGG analysis was used to explore the biological pathways.

Based on the clinical information of PAAD patients downloaded from TCGA database, the survival time and state were extracted. Kaplan–Meier survival analysis and univariate Cox regression analysis were conducted by integrating clinical data and ceRNA network node information to evaluate the prognostic value of node genes. Lasso regression (conducted via the “glmnet” package in R language) was used to remove redundancy factors to prevent over-fitting of the model, and then the ceRNA network prognosis model was obtained through multivariate Cox regression analysis.

All patients were divided into high-risk and low-risk groups according to the median of the model’s risk score, and the survival curve was drawn to compare the differences in survival between them. The area under the curve (AUC) of the receiver-operating characteristic (ROC) curve (visualized by the “plotAUCcerve” package in R language) was used to evaluate the prediction efficiency of the model. Next, a decision curve analysis (DCA) curve was plotted with the “ggDCA” package of R, and a nomogram was constructed with the “rms” package of R. Furthermore, the calibration curve was used to evaluate the accuracy of the nomogram in predicting the survival rate of patients with PAAD.

The Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org) (Yang et al., 2013), the largest available database, was used to evaluate the sample’s half-maximal inhibitory concentration (IC50) and predict the chemotherapy response of samples with the “pRRophetic” package of R. We selected eight common medicines (A.443654, BI.2536, lapatinib, dasatinib, gemcitabine, docetaxel, paclitaxel, and cyclopamine) to predict the IC50 value of each sample. Subsequently, the difference between the IC50 values in the high-risk and low-risk groups was analyzed using the Wilcoxon rank-sum test. Additionally, the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm and subclass mapping were used to predict the immune responses of the sample to anti-PD-1 and anti-CTLA-4 treatment.

The abundance of 22 immune cell subgroups in 178 tumor samples was calculated by the “CIBERSORT” package in R language, with p < 0.05 as the threshold. The “Pheatmap” package of R was used to visualize the infiltration of immune cells in the normal and tumor groups. The Wilcoxon rank-sum test was used to analyze the differentially expressed tumor-infiltrating immune cells in PAAD tissues and normal tissues. Then, Kaplan–Meier analysis was conducted to identify prognosis-related immune cells.

We used Kaplan–Meier survival analysis and univariate Cox regression to investigate the survival-related immune cells. At the same time, the final immune cell model was obtained by Lasso regression and multivariate Cox regression analysis, and the model’s accuracy was evaluated by the AUC value. We used the multi-index AUC value and DCA to estimate the quantification of the prognostic model and utilized the calibration curve to assess the predictive accuracy of the nomogram. Meanwhile, we used Pearson correlation analysis to explore the relationship between key genes in the ceRNA network and immune cells in the model.

Hereby, we confirmed the tissue and cell expression levels of one lncRNA and three mRNAs in the model through experiments. In particular, TCGA and GTEx RNAseq data in TMP format that were extracted from the UCSC Xena browser (https://xenabrowser.net/datapages/) (Vivian et al., 2017) were uniformly processed using the Toil process. The PAAD data in TCGA and the corresponding normal tissue sample data in GTEx were extracted for differential expression analysis with the “EdgeR” package of R. Then, the expression levels of key genes were validated by real-time quantitative reverse transcription PCR (RT-qPCR) in PAAD tissues and cell lines.

To extract total RNA from tissue, soybean-sized samples of PAAD tissues and paracancerous tissues were removed from liquid nitrogen and then immediately treated with Trizol (Invitrogen, Carlsbad, CA, United States). Then, RT-qPCR was applied using the cDNA Kit (Genstar, Beijing, China) and qPCR Kit (Genstar, Beijing, China) based on the manufacturer’s instruction. The relative mRNA expression was calculated with the 2^–ΔΔCT method. The primers are presented in Table 1.

All statistical analyses were performed with R (version 4.0.2). Bioconductor packages (Eivazi et al., 2016) used were as follows: limma, GDCRNATools, edgeR, ggplot2, survminer, timeROC, ggDCA, rms, glmnet, plotAUCcerve, and preprocessCore. Co-expression analysis was performed using the Pearson correlation analysis. The comparison of means was conducted with the independent sample Student’s t-test. GraphPad Prism 8 was used for the corresponding statistical analysis. p values < 0.05 on two sides were considered statistically significant. The flow chart of this study is shown in Figure 1.

A total of 412 DEGs in PAAD were identified from TCGA database with the filter threshold of FC > 2 and false discovery rate (FDR) < 0.01, including 360 differentially expressed mRNAs (147 downregulated and 213 upregulated; Figures 2A,B), 28 differentially expressed miRNAs (16 upregulated and 12 downregulated; Figures 2C,D), and 24 differentially expressed lncRNAs (14 upregulated and 10 downregulated; Figures 2E,F).

After miRcode matching and miRNA target gene prediction, 12 differentially expressed lncRNAs, 13 miRNAs, and 145 mRNAs were identified. With 132 nodes and 311 edges, the ceRNA regulatory network of the PAAD-related DEGs for lncRNA, miRNA, and mRNA was constructed (Figure 3A). Function analysis of the KEGG signaling pathway and GO analysis (Figures 3B,C) revealed that the DEGs were mainly enriched in the integrin-mediated signaling pathway and chemokine signaling pathway, with a close relationship to regulated exocytosis, actin cytoskeleton organization, regulation of response to biotic stimulus, and regulation of leukocyte-mediated immunity. Furthermore, 11 significantly correlated GO and KEGG items were determined and visualized according to genes (CCs, MFs, BPs, and KEGG signaling pathway; Figure 3D).

Kaplan–Meier survival analysis and log-rank tests were used to determine the survival difference of genes in the ceRNA network, which identified 41 genes correlated with PAAD prognosis. Then, the top 12 genes were filtered according to their p value with the log-rank test, and the Kaplan–Meier survival curve was plotted (Figure 4). Key genes in the network were further identified by univariate Cox regression analysis and Lasso regression analysis. We found that nine significant genes were included in the multivariate Cox regression analysis (Figures 5A,B). Next, through multivariate Cox regression analysis, five key genes in the ceRNA network were used to establish the prognostic model of PAAD, namely, LRRC1, RNF166, LY6D, MIR600HG, and hsa-miR-424-5p (Figure 5D and Table 3). The results of Kaplan–Meier curves indicated that patients in the high-risk group obtained a poorer prognosis than those in the low-risk group (Figure 5C). Meanwhile, the AUC of the ROC curve was generated to assess the predictive accuracy of the Cox regression risk model. Next, we found that its accuracy was great (the AUC values for one-, three-, and five-year survival were 0.741, 0.819, and 0.804, respectively) (Figure 5E). The comparison of multi-index AUC values at multiple time points revealed that the prediction accuracy of the early-stage model was better than that of the individual genes, with hsa-miR-424-5p at the end stage showing better prediction (Figure 5F). The DCA indicated that the risk model had better prediction performance than the model constructed by traditional clinicopathological features (Figure 5G). The final nomogram was constructed by the expression of key genes in the model, and the calibration curve illustrated the higher predictive accuracy of the nomogram (Figures 5H,I).

The chemotherapeutics gemcitabine (Figure 6A), dasatinib (Figure 6B), lapatinib (Figure 6C), docetaxel (Figure 6D), cyclopamine (Figure 6E), paclitaxel (Figure 6F), A.443654 (Figure 6G), and BI.2536 (Figure 6H) were compared in the high-risk and low-risk groups. We found that gemcitabine had no difference between the high-risk and low-risk groups, and the remaining seven drugs showed higher sensitivity in the high-risk group than the low-risk group. As shown in Figure 6I, low-risk patients had a significantly positive response to anti-PD-1 treatment. Differential analysis of immune checkpoint expression revealed that the expression levels of immune checkpoint factors such as PD-1, CTLA-4, LAG3, and BTLA were significantly different between the high-risk and low-risk groups (Figure 6J).

The percentage of 22 immune cells in all samples was obtained using the deconvolution algorithm CIBERSORT, and a histogram and heatmap of the immune cells were drawn (Figures 7A,B). We found that naive CD4 T cells, resting memory CD4 T cells, eosinophils, and gamma delta T cells were significantly infiltrated in PAAD samples, whereas memory B cells, memory monocyte B cells, and resting dendritic cells were significantly infiltrated in normal tissue samples (Figure 7C). Moreover, Pearson correlation analysis suggested that naive CD4 T cells were significantly positively correlated with memory B cells (r = 0.58), and activated NK cells were significantly negatively correlated with resting NK cells (r = −0.58) (Figure 8A). The PAAD prognosis-related immune cells were identified by Kaplan–Meier survival analysis and log-rank tests, which showed that the infiltration levels of memory B cells and activated dendritic cells in tumor tissues were significantly correlated with PAAD prognosis (Figures 8B,C). Subsequently, under Lasso regression (Figures 9A,B) and Cox regression analysis (Figure 9D, Table 2), CD8 T cells, M1 macrophages, and activated dendritic cells were incorporated into the model. The results of Kaplan–Meier curves demonstrated that patients in the high-risk group obtained a poorer prognosis than those in the low-risk group (Figure 9C). The levels of CD8 T cell, M1 macrophage, and activated dendritic cell infiltration in the high-risk and low-risk groups are shown via a heatmap (Figure 9F). The AUC of the ROC curve was generated to assess the predictive accuracy of the Cox regression risk model, which was found to be excellent (the AUC values for one-, three-, and five-year survival were 0.766, 0.768, and 0.669, respectively) (Figure 9E). The comparison of the multi-index AUC values at multiple time points revealed that the predictive efficacy of the model was significantly better than that of individual immune cells at all time points (Figure 9G). The DCA indicated that the immune risk model had better predictive performance than the model constructed by traditional clinicopathological features (Figure 9H). The final nomogram was constructed based on the content of immune cells in the model, and the calibration curve illustrated the high accuracy of the nomogram in predicting survival (Figures 9I,J).

Co-expression analysis was conducted through the Pearson correlation test to explore the correlation between key genes in the ceRNA network and prognosis-related immune cells (Figure 10A). We found that MIR600HG and M1 macrophages were inversely correlated (R = −0.250; p < 0.01) (Figure 10B), MIR600HG and CD8 T cells were positively correlated (R = 0.270; p < 0.01) (Figure 10C), RNF166 and CD8 T cells were positively correlated (R = 0.280; p < 0.01) (Figure 10D), and LRRC1 and M1 macrophages were significantly positively correlated (R = 0.355; p < 0.001) (Figure 10E).

Compared to normal tissues, we found that the expression levels of LRRC1, RNF166, and LY6D were significantly higher in PAAD tissues (Figures 11A–C), while the expression of MIR600HG was significantly lower in PAAD tissues (Figure 11D). Meanwhile, RT-qPCR revealed that the expression levels of LRRC1, RNF166, and LY6D were increased in PAAD tissues; however, the expression of MIR600HG was decreased in PAAD tissues (Figure 11I). Excitingly, similar results were observed in PAAD cell lines compared to HPDE cell lines (Figures 11E–H).