Transcriptome expression profiles, mutation data, and relevant clinical information of patients with PAAD were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository). The data were then collated using Strawberry Perl (version 5.32.1.1) scripts to obtain mRNA and lncRNA data matrixes for subsequent studies. The 86 CAFs related genes used for the study were obtained from The Human Gene Database (https://www.genecards.org/), with a relevance score of >11 (Supplementary Appendix S1).

The mRNA expression data of 86 CAFs-related genes were extracted using the R package “limma.” The set of CAFs-related lncRNAs was obtained by co-expression analysis using a Pearson correlation coefficient >0.4 and a threshold of p < 0.001. The correlation data of CAFs-related genes and lncRNAs were constructed using the R packages “dplyr,” “ggalluvial,” and “ggplot2,” and correlation Sankey plots were generated. Subsequently, the “limma” package was used to merge the survival and lncRNA expression data of each PAAD patient.

The R packages “caret,” “timeROC,” “survminer,” “survival,” and “glmnet”were used to establish the risk signature of CAFs-related LncRNAs in PAAD. Data of the patients obtained from the TCGA were randomly divided into training and validation groups by 1:1, and various clinical traits were analyzed in the two groups. Then the prognosis-associated lncRNAs in the training group were obtained using univariate Cox analysis, and the prognostic forest plots (p < 0.05) were plotted. Subsequently, the “pheatmap” package was used to plot the expression heat map of prognosis-related lncRNAs. The LASSO regression analysis was used to screen candidate lncRNAs to avoid overfitting. Risk scores for all patients were obtained by the following formula: risk score = (coefficient lncRNA1*lncRNA1 expression) + (coefficient lncRNA2*lncRNA2 expression) + …+ (coefficient lncRNAn* lncRNAn expression). Coefficient and expression represent the regression coefficient and expression values of the CAFs-related lncRNAs model, respectively. Based on the median risk score in the training cohort, patients in the validation and training cohorts were split into high- and low-risk groups. In addition, the R packages “tidyverse,” “ggplot2,” and “ggExtra” were used for expression correlation analysis of CAFs-related genes and model lncRNAs and to draw a correlation heat map. The “limma,” “survivor,” and “survminer” packages were used to plot Kaplan-Meier (KM) curves for the high- and low-expression groups of model lncRNAs.

The R packages “pheatmap,” “survival,” and “survminer” were used for survival analysis of the training, validation, and entire cohorts. Risk curves, risk heat maps, and survival status maps were plotted for each cohort. In addition, overall survival (OS) and progression-free survival (PFS) survival curves were also plotted for the all cohort. Univariate and multivariate Cox regression analyses were used to assess whether the risk score and selected clinical characteristics were independent prognostic factors. Additionally, ROC curve analyses were performed using “survminer”, “survival”, and “timeROC” to assess the prognostic value of the developed signature by the area under the curve. Furthermore, the “survival,” “rms,” and “pec” packages were used to calculate the consistency index (c-index) to evaluate the best prediction of the model.

The R package “ComplexHeatmap” was used to plot the heat map of the relationship between the high- and low-risk groups of the model and different clinicopathological parameters. Meanwhile, “survminer” and “survival” were used to plot the survival curves of high- and low-risk groups to determine whether the constructed risk model applied to PAAD patients with different clinicopathological parameters.

Based on the results of the multivariate regression analysis, risk status and age were used to construct nomograms for 1-year, 3-years, and 5-years OS using the “regplot,” “survival,” and “rms” R packages. Hosmer-Lemeshow test calibration curve (method = “boot”, B = 1,000) was utilized to validate if the actual results correlate with the anticipated results.

The Gene Set Variation Analysis (GSVA) computationally detects differences in pathway activity in a sample population (Hänzelmann et al., 2013). GSEA enrichment analysis by the R package “limma,” “GSEABase,” and “GSVA” to obtain the enrichment of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway in the high-risk and low-risk groups. The R package “pheatmap” was used to plot the pathway enrichment heatmap. In addition, differentially expressed genes between high- and low-risk populations (A fold change >2 and FDR <0.05) were identified using the R package “limma.” Subsequently, “enrichplot,” “GOplot,” “ggplot2,” “org.Hs.eg.db,” and “clusterProfiler” packages were used to complete GO analysis and explore potential pathways.

The Strawberry Perl script collated the PAAD mutation data downloaded from the TCGA to obtain the tumor mutation burden (TMB) data for each patient. R packages “Limma” and “ggpubr” were used to analyze the TMB differences between the high- and low-risk groups of the model and generate violin plots. Meanwhile, “maftools” was used to map the mutation waterfall of the 15 genes with the highest mutation frequency in the high- and low-risk groups. In addition, the R software was used to obtain the optimal cutoff values of tumor mutation burden and classify patients into low-TMB and high-TMB groups. The “survivor” and “survminer” packages were used to plot the survival curves of patients in two risk groups and the survival curves of patients in the high- and low-TMB groups combined with the high- and low-risk groups.

CAFs play an essential regulatory role in the tumor immune microenvironment. To further explore the correlation between risk models constructed by CAFs-related lncRNA signature and immune microenvironment, the R packages “ggtext,” “tidyverse,” “ggpubr,” “scales,” and “ggplot2” were used to analyze the correlation between immune cells and risk scores and generate correlation bubble plots. In addition, single sample GSEA (ssGSEA) was performed to classify gene sets with common physiological regulation and biological functions (Subramanian et al., 2005). The “GSEABase” and “GSVA” packages were utilized for ssGSEA analysis to calculate the immune cell and immune-related function scores of the samples. “reshape2,” “ggpubr”, “pheatmap,” and “reshape2” packages were used to obtain box plots for ssGSEA differential analysis and heat maps for immune function differences in high- and low-risk groups of the model.

The potential predictive value of CAFs-related lncRNA signature for immune checkpoint efficacy was explored by utilizing the “ggplot2” and “ggpubr” packages to analyze the differences of immune checkpoint-related genes between the high- and low-risk groups and generate box plots of differentially expressed genes. Immune checkpoint programmed death-ligand 1 (PD-L1) on cancer cells binds to programmed cell death-1 (PD-1) on immune cells and contributes to the immune escape of tumor cells (Yi et al., 2021). Furthermore, in some cases, tumor PD-L1 expression correlates with immunotherapy response (Yu et al., 2016; Yi et al., 2018). “ggpubr” and “limma” packages analyzed the correlation between PD-L1, PD1, and CTLA4 expression and seven model lncRNAs and “corrplot” package plotted the correlation.

ESTIMATE is a new algorithm for counting immune and stromal cells infiltrating tumor tissue (Yoshihara et al., 2013). Using the R packages “ESTIMATE” and “limma,” the amount of immune and stromal cells in the tumor tissue of each PAAD case was evaluated to determine the corresponding scores. The sum of the stromal and immune scores is the ESTIMATE score, which is inversely linked to tumor purity. Box plots illustrating the differences in stromal, immune, and ESTIMATE scores between high- and low-risk groups were generated using the “ggpubr” package.

The potential clinical significance of CAFs-related lncRNA signatures in chemotherapy and targeted therapies were explored using the half-maximal inhibitory concentrations (IC50) of different drugs in the high-risk and low-risk groups, obtained using the R packages “ggpubr” and “pRRophetic” (Geeleher et al., 2014). Drugs with different IC50s in the two groups were represented as box plots (p < 0.001).

We obtained expression data from TCGA database for 179 PAAD tumor samples. By co-expression analysis of mRNAs of CAFs-related genes, we obtained 378 CAFs-associated lncRNAs (correlation coefficient >0.4, p < 0.001) (Figure 1A).

The TCGA cohort was divided into training and validation groups, and clinicopathological characteristics were compared (Table 1). Univariate Cox (uni-Cox) regression analysis obtained 72 CAFs-related lncRNAs in the training group that was significantly associated with the OS of patients (p < 0.05). Figures 1B,C show the prognostic forest plot and expression heat map of the 72 LncRNAs, respectively. We further performed lasso regression analysis (Figures 2A,B) and extracted seven of these LncRNAs for model construction (Table 2). Risk scores were calculated from lncRNAs screened by lasso regression to create the formula: Risk score = AP005233.2 × (0.257) + AC090114.2 × (−0.829) + DCST1-AS1 × (0.596) + AC092171.5 × (−0.561) + AC002401.4 × (0.215) + AC025048.4 × (−1.388) + CASC8 × (0.348). The correlation between the expression of seven lncRNAs and CAFs-related genes was demonstrated with a heat map (Figure 2C).

Survival analysis showed that patients with high expression of AC090114.2, AC092171.5, and AC025048.4 had significantly better survival rates than the low expression group, while the remaining four genes showed the opposite effect (Figures 2D–J). The expression heat map showed that AC090114.2, AC092171.5, and AC025048.4 were low in both the training, validation, and entire cohorts in the high-risk group. Whereas consistent with the above results, the remaining four genes showed the opposite effect (Figures 3A–C). Subsequently, the survival status, risk score distribution, and OS survival curves of the low- and high-risk patients in the training, validation, and entire cohorts were assessed using the risk scores. All results showed that the prognosis of the low-risk group was significantly better than that of the high-risk group (Figures 3D–L). Similar results were obtained for the PFS survival curves of the training and entire TCGA cohorts. Although there was no statistical difference in PFS between the high- and low-risk groups in the validation group, a trend towards a separation of survival curves was observed (Figures 3M–O).

A heat map of clinicopathological parameters showed differences in tumor grade between high- and low-risk groups (Figure 4A). Survival analysis showed that PAAD patients with different gender, ages, tumor grades, and stages all survived significantly better in the low-risk group than in the high-risk group (Figures 4B–I), demonstrating the model’s applicability to PAAD patients with different clinicopathological parameters. Further, univariate Cox (uni-Cox) regression and multivariate Cox (multi-Cox) regression suggested the risk score as an independent prognostic factor with hazard ratios (HR) of 1.174 and 1.190, respectively, with 95% confidence intervals (CI) of 1.124–1.227 (p < 0.001) and 1.135–1.248 (p < 0.001) (Figures 5A,B). Also, the patient’s age was an independent prognostic parameter. In addition, the ROC curve was used to assess the sensitivity and specificity of the risk model to the prognosis of PAAD. The results showed that the model had a significantly higher predictive value than other clinicopathological parameters, with an area under the curve (AUC) of 0.811, 0.816, and 0.840 at 1, 3, and 5 years respectively (Figures 5C–F). Furthermore, the c-index of the risk score was also higher than that of the other clinical parameters (Figure 5G). Together, these results demonstrate the good performance of the risk model.

We constructed a nomogram based on the patient’s age and risk status to facilitate prognosis prediction for PAAD patients (Figure 5H). The corresponding scores for the patient’s age and risk status were calculated in the nomogram, and the total score was used as a prognostic prediction tool. A calibration curve was also plotted (Figure 5I). The results showed a good agreement between the survival of the PAAD patients and the values predicted by the nomogram.

To explore the differences in biological behavior between high- and low-risk groups, we used GSVA to investigate the differences in functional pathways between the groups. The pathways enriched in the high-risk group included the cell cycle, DNA replication, p53 signaling pathway, mismatch repair, and fatty acid metabolism, which were associated with tumor invasion. On the other hand, functions such as intestinal immune network, chemokine signaling pathway, and glycan degradation were enriched in the low-risk group (Figure 6A). We further analyzed the enrichment of differentially expressed genes (DEGs) in different risk groups in terms of biological functions by GO analysis. The results suggested that DEGs were enriched in functions such as cellular ion channels, membrane receptors, T-cell receptors, and transduction of signals (Figures 6B,C).

Tumor mutation burden (TMB) is defined as the number of somatic mutations per megabase. TMB is a crucial driver in generating immunogenic neopeptides in tumor cells and affecting the patient’s response to ICIs (Sha et al., 2020). We used the TGCA somatic mutation data to generate TMB scores. Further analysis revealed that TMB levels were significantly lower in the low-risk group than in the high-risk group (Figure 7A). Also, survival analysis showed that higher TMB in PAAD was associated with a poorer OS (Figure 7B). Given the prognostic role of the risk model and TMB in PAAD, we further explored the prognostic value of combining the two by dividing all samples into four groups: high-TMB/high-risk, low-TMB/low-risk, high-TMB/low-risk, and low-TMB/high-risk. The results showed a significant difference in survival between the four groups (p < 0.001), with patients with high-TMB/high-risk having the worst OS and those in the low-TMB/low-risk group having the best overall survival (Figure 7C). In addition, the frequency of mutations was higher in the high-risk group (96.88%) than in the low-risk group (62.12%) (Figures 7D,E). The highest mutation frequencies were found in KRAS (79%), TP53 (68%) and SMAD4 (24%) in the high-risk group, while the highest mutation frequencies in the low-risk group were in TP53 (39%), KRAS (35%) and SMAD4 (18%).

The tumor immune microenvironment (TIME) is closely related to the prognosis of patients with tumors. There is growing evidence that CAFs can synergize with TIME components, particularly immune cells, to render the tumor microenvironment (TME) immunosuppressive (Barrett R. and Puré E., 2020; Barrett R. L. and Puré E., 2020). To further explore the correlation between the CAFs-related lncRNAs signature and TIME, we analyzed the association between immune cells and risk scores. The results indicated that most immune cells were negatively related with the risk scores. In contrast, regulatory T cells (Tregs) on the CIBERSORT platform and common lymphoid progenitor, CD4⁺ Th1/2 cells on the XCELL platform, were positively correlated with the risk scores (p < 0.05) (Figure 8A). In addition, ssGSEA analysis showed that CD8⁺ T cells, dendritic cells (DCs), plasmacytoid dendritic cells (pDCs), neutrophils, mast cells, helper T cells, and tumor-infiltrating lymphocytes (TIL) were significantly fewer in the high-risk group than in the low-risk group (Figure 8B). In terms of immune-related functions, T cell co-stimulation, T cell co-inhibition, type II interferon (IFN) responses, and cytolytic activity were significantly weaker in the high-risk group than in the low-risk. The opposite trend was seen with the major histocompatibility complex (MHC) class I and type I IFN response (Figures 8C,D). Further, we analyzed the relationship between risk groups and ESTIMATE scores. The results suggested that the ESTIMATE score, stromal score, and immune score were significantly higher in the low-risk group than in the high-risk group (Figures 8E–G), which corroborated with the above results and together indicated that the low-risk group had a higher immune cell infiltration status.

ICIs therapy offers a new tool for clinical cancer treatment by enhancing anti-tumor immune responses through a regulatory pathway of T cells (Sharma and Allison, 2015). However, immune checkpoint therapy benefits only a small proportion of patients with specific tumor types, and one of the main problems is the lack of validated prognostic biomarkers (Sharma et al., 2021). The current ICIs target the programmed death-1 (PD-1), programmed death-ligand 1 (PD-L1), and cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) ((Sharma et al., 2021). Although the expression of PD-L1 did not differ in the high- and low-risk groups, the mRNAs of most other immune checkpoint-related genes, including PD1 and CTLA4, were highly expressed in the low-risk group (Figure 8H). Our results suggest that low-risk patients may benefit more from ICIs. Finally, we analyzed the correlation of the expression between the seven lncRNAs and the three immune checkpoints. The results showed that DCST1-AS and AC092171.5 were negatively correlated with the expression of all three immune checkpoints (p < 0.05, Figures 8I–K). Together, these results suggest that CAFs-related lncRNAs signature might better distinguish PAAD patients with different tumor immune microenvironment characteristics and provide a basis for selecting clinical immunotherapy.

Using the pRRophetic algorithm analysis, we found that the IC50 of some compounds differed between the high- and low-risk groups of the model (p < 0.001) (Figures 9A–P). Among them, the IC50 of mTOR inhibitor AZD8055, All-trans retinoic acid (ATRA), lestaurtinib (CEP-701), Bcl-2 family protein inhibitor navitoclax (ABT-263), and other drugs were higher in the high-risk group than in the low-risk group. The opposite effect was shown by drugs such as doxorubicin, gefitinib, and mitomycin C.