The public RNA sequencing of 182 patients and full clinical information were obtained from the GDC data portal using gdc-client (https://gdc-portal.nci.nih.gov/). The clinical characteristics including (age, gender, clinical grade and clinical stage) and survival overall survival (OS) information. Patients without survival information were excluded from further evaluation. In addition, the expression level of 60 ferroptosis-related genes (ACSL4, AKR1C1, AKR1C2, AKR1C3, ALOX15, ALOX5, ALOX12, ATP5MC3, CARS1, CBS, CD44, CHAC1, CISD1, CS, DPP4, FANCD2, GCLC, GCLM, GLS2, GPX4, GSS, HMGCR, HSPB1, CRYAB, LPCAT3, MT1G, NCOA4, PTGS2, RPL8, SAT1, SLC7A11, FDFT1, TFRC, TP53, EMC2, AIFM2, PHKG2, HSBP1, ACO1, FTH1, STEAP3, NFS1, ACSL3, ACACA, PEBP1, ZEB1, SQLE, FADS2, NFE2L2, KEAP1, NQO1, NOX1, ABCC1, SLC1A5, GOT1, G6PD, PGD, IREB2, HMOX1, ACSF2) were constructed from the dataset.

Based on the lncRNA annotation file downloaded from the GENCODE (https://www.gencodegenes.org/human/) website and Ensemble IDs, we annotated 14,086 lncRNAs and 19,604 mRNAs according in the TCGA datdset. Pearson correlation analysis was used to evaluate the ferroptosis-related lncRNA. LncRNAs with p value <0.01 and an absolute Pearson correlation coeffcient≥ 0.4 or ≤0.4 were selected as ferroptosis-related lncRNAs. Then, the univariate cox regression analysis was used to screen out the prognostic lncRNAs(p < 0.05).

The risk model was constructed through the Least Absolute Shrinkage and Selection operator (LASSO) and multivariate Cox regression analysis. To improve the accuracy of the statistical model, Lasso Cox regression analysis was carried out by using the R package “glmnet”. There are 10 ferroptosis-related lncRNAs were screened to construct the best prognostic model through the multivariate Cox regression analysis. The risk score was generated using the following formula: risk score  =   ∑ CoeflncRNAs   × ExplncRNAs All patients were divided into high-risk and low-risk groups according to the median cutoff of risk score.

The Kaplan-Meier survival method was applied to evaluate the availability of the prognostic model and all the statistical analyses were conducted using R language (version 4.0). In addition, the area under the ROC curve (AUC) was calculated to evaluate the prognostic accuracy and sensitivity of the model. Survival analysis was also performed for the five crucial lncRNAs(ZNF236-DT, CASC8, PAN3-AS1, SH3PXD2A-AS1, LINP1) respectively in the model. We divided pancreatic cancer patients into different subgroups and compared whether there were differences in risk scores among the groups.

Nomogram is a tool that can personally calculate the survival rate of patients with specific tumors and has great value in clinical application (Park, 2018). So, we constructed a nomogram to help clinicians conveniently use our model to predict the overall survival at 1, 3, and 5 years of patients with pancreatic cancer. The nomogram is completed by use the R package “rms”, “Hmisc”, it includes risk score, age, gender, grade and stage.

We identified differentially expressed genes (DEGs) between low-risk and high-risk subgroups in the TCGA cohort in order to investigate potential biological functions and pathways between the two groups. All the DEGs in the study were performed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis by using the clusterProfile R package, which the results with p < 0.05 were considered significant.

The proportion of tumor-infiltrating immune cells between two risk groups were calculated by CIBERSORT, ESTIMATE, single-sample gene set enrichment analysis (ssGSEA) algorithms based on FRlncRNAs signature. The results were filtered with p value < 0.05.

Immune checkpoints are the immunosuppressive pathways that immune cells possess to regulate and control the persistence of the immune responses (Abril-Rodriguez and Ribas, 2017). The up-regulation of immune checkpoint protein is also one of the mechanisms by which tumor cells evade immune responses (Pardoll, 2012). In recent years, checkpoint blockade therapies have become a new type of cancer treatment. Thus, we investigated whether there were differences in the expression of immune checkpoint genes between the two risk groups, which aimed to investigate the potential role of FRlncRNAs signature in immune checkpoint blockade therapy.

Finally, we used a variety of currently acknowledged methods to analyze the relationship between risk and immune-cell characteristics, including XCELL, TIMER, QUANTISEQ, MCPOUNTER, EPIC, CIBERSORT-ABS and CIBERSORT. The relationship between the risk score and the immune infiltrated cells was performed by Spearmen correlation analysis. The correlation coefficient is bounded by 0 and the results were shown in a lollipop diagram.

Human pancreatic cancer cells PANC-1 were used in our experiment. Cells were cultured in high-glucose Dulbecco’s Modified Eagle’s Medium (supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin) at 37°C and 5% CO2 in a cell incubator. The siRNA against SH3PXD2A-AS1 (siSH3PXD2A-AS1), the LINP1(siLINP1) and the control siRNA (siNC) were purchased from RiboBio Corporation (Guangzhou, China). The siRNAs transfection were performed using Lipofectamine 2000 transfection reagent. The sequence of siLINP1 and siSH3PXD2A-AS1 were showed below:

siLINP1

5′-CCA​ACT​GCG​GGA​CTT​CAG​A dTdT-3′(sense), 5′-TCT​GAA​GTC​CCG​CAG​TTG​G dTdT- 3′(antisense).

siSH3PXD2A-AS1

5′-CCC​TAA​GGA​CAG​AAT​GCA​A dTdT-3′(sense), 5′- TTG​CAT​TCT​GTC​CTT​AGG​G dTdT- 3′(antisense).

PANC-1 cells were grown in six-well plates until they formed a tight cell monolayer. A 200ul sterile tip was used to form a wound in the middle of the cell. Cells were then cultured in medium supplemented with 2% FBS. The width of wounds were photographed using a microscope at 0 and 24 h. The percentage of wound healing was generated using the following formula: average of ([gap area: 0 h]–[gap area: 24 h])/[gap area: 0 h]).

Firstly, PANC-1 cells in 200 µl of serum-free DMEM were incubated to the upper chamber which is coated with Matrigel, 500 µl of complete medium was placed in the lower chamber. The upper chambers containing non-invasive cells were removed with cottons swabs after 48 h. Then, the lower chamber was fixed with 4% paraformaldehyde and stained in crystal violet. The chambers were observed under microscope.

For Cell Counting Kit-8 (CCK-8) assays, cells transfected for 48 h were collected and seeded into a 96-well plate. According to the instructions, the optical density (OD) values were measured at a 450-nm wavelength after the cells were incubated for 2 h at 37°C.

There are a total of 182 patients from the TCGA dataset who were identified and included in this study. 14,086 lncRNAs in the TCGA datasets were obtained by matching the lnRNA annotation with the ENSEMBL ID downloaded from the “GENCODE” website. 60 ferroptosis-related genes were collected from the published literature for subsequent analysis. LncRNAs with p value < 0.01 and an absolute Pearson correlation coeffcient ≥ 0.4 or ≤ 0.4 were selected as ferroptosis-related lncRNAs. Among 246 ferroptosis-related lncRNAs, 89 lncRNAs were identified as prognostic ferroptosis-related lncRNAs through the univariate cox regression analysis (p < 0.05). The research process is shown in Figure 1A and the forest map of these 89 lncRNAs is shown in Figure 1B.

To build the FRLS for the patients of pancreatic cancer in TCGA dataset, the above 89 prognostic ferroptosis-related lncRNAs were incorporated into the least absolute shrinkage and selection operator (LASSO) regression firstly (Figures 1C,D). 10 ferroptosis-related lncRNAs including ZNF236-DT, AC005332.6, AC019186.1, AC087501.4, CASC8, PAN3-AS1, LINC01133, SH3PXD2A-AS1, LINP1, AC090114.2 were generated. To further optimize the results, we performed the multivariate Cox regression analysis to construct a prognostic model for OS using the expression level of the 10 ferroptosis-related lncRNAs. An optimal 5 lncRNAs (ZNF236-DT, CASC8, PAN3-AS1, SH3PXD2A-AS1, LINP1) signature and coefficient of each were identified (Figures 1E,F). The forest plot shows that ZNF236-DT and PAN3-AS1 are protective factors with HR (Hazard Ratio) < 1, while CASC8, SH3PXD2A-AS1 and LINP1 are risk factors with HR > 1 (Figure 1E). The patients were further divided into high-risk group (n = 88) and low-risk group (n = 89). The groups were generated by the formula mentioned in section “Materials and Methods”. Principal Component Analyses (PCA) showed a significant distribution difference between high- and low-risk subgroups (Figures 1G,H).

The Kaplan-Meier survival analysis showed that the high-risk group (n = 88) had a lower probability of survival and a shorter overall survival (OS) time than patients in the low-risk group (n = 89) (Figure 2A). The area under ROC curve (AUC) reached 0.733, which demonstrated that the FRLS could be used to predict OS in the TCGA cohort (Figure 2C). Risk score, survival status and heatmap are plotted in Figure 2B. Then, the 5 crucial lncRNAs(ZNF236-DT, CASC8, PAN3-AS1, SH3PXD2A-AS1, LINP1) of the FRLS were evaluated by survival analysis. The Kaplan–Meier survival curves confirmed the forest plot shown in Figure 1E, which means higher expressions of ZNF236-DT, PAN3-AS1 and lower expressions of CASC8, SH3PXD2A-AS1 and LINP1 were associated with better overall survival (OS) (Figures 2D–H).

To further assess the prognostic efficacy of FRLS, we divided patients into various subgroups based on clinicopathological features and compared the levels of risk score between different groups. No risk score differences were observed between patients satisfied by gender and age (age was divided by 65 years old) (Figures 3A,B). But the results suggested that patients with the clinicopathological characteristics of grade 2 and 3, stage II have higher levels of risk score compared with other corresponding subgroups (Figures 3C,D). Thus, FRLS may have certain predictive value for clinicopathological features. The heatmap in Figure 3E further shows the distribution of clinicopathological characteristics and the risk groups.

Afterward, we have used univariate and multivariate Cox analyses to identify independent prognostic factors indicator. The results of univariate Cox and multivariate Cox regression analysis showed that only the risk score was significantly associated with OS in TCGA dataset (univariate Cox analyses: HR = 1.412, 95% CI = 1.286–1.551, p < 0.05) (Figure 3G). As shown in Figure 3F, the AUC value of the risk score was 0.724, higher than the AUC values of other clinicopathological factors. These results confirmed the FRLS is an independent and reliable prognostic indicators for pancreatic cancer.

In summary, we constructed a nomogram using the risk score (based on FRLS) and clinicopathological features, including age, gender, WHO grade and stage in the TCGA dataset (Figure 3H).

To investigate potential biological functions and pathways between the two risk subgroups, we performed functional enrichment analysis of the differentially expressed genes (DEGs) between the two groups. There are 973 DEGs [|log2 (fold change)| > 1 and p < 0.05] between the low- and high-risk subgroups. Interestingly, the GO analysis revealed that the DEGs were highly enriched in several immune-related biological processes, including immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, humoral immune response, lymphocyte mediated immunity, complement activation, B cell receptor signaling pathway (Figure 4A). Similarly, the KEGG analysis also showed significant enrichment of DEGs in immune-related pathways, for instance, cytokine-cytokine receptor interaction, chemokine signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway (Figure 4B).

Since the KEGG and GO analysis suggested that DEGs were enriched in immune-related functions and pathways, we further investigated the immune landscape in TCGA dataset. The ESTIMATE algorithm showed that the immune score, stromal score and ESTIMATE score were all higher in the low-risk group than in the high-risk group (Figures 5A–C). These results also suggest that tumor purity was lower in the low-risk group than in the high-risk group.

Moreover, the CIBERSORT algorithm was used to analyze the 22 different immune cell types among the two risk group. As shown in Figures 5B,D cell naive, T cell CD8 and T cells regulatory (Tregs) were up-regulated in the low risk subgroup of the TCGA cohort, while NK cells activated and mast cell resting were significantly down-regulated (p < 0.05). Next,we further analyzed the correlation between FRLS and immune cells. The results indicated that this signature was most significantly positive correlation with immune infiltration of Macrophages M0(Cor = 0.25, p = 0.004), Macrophages M1(Cor = 0.2, p = 0.022), Macrophages M2(Cor = 0.22, p = 0.012) and Mast cells resting (Cor = 0.2, p = 0.02), but negatively correlated with B cells naive (Cor = −0.35, p = 4.5e-05), Plasma cells (Cor = −0.2, p = 0.023), T cells CD4 memory activated (Cor = −0.2, p = 0.022), T cells CD8 (Cor = −0.3, p = 0.00059) and T cells regulatory (Tregs, Cor = −0.27, p = 0.0014) (Figures 5F–N). These findings again confirmed that this FRlncRNAs signature was associated to immune cell infiltration in pancreatic cancer.

By using the ssGSEA analysis, we quantified the enrichment scores of immune-related functions in the two risk groups. The results showed that the high-risk group was significantly correlated with most immune-related functions, such as CCR, Check-point, Cytolytic-activity, HLA, Inflammation-promoting, T cell co-inhibition, T cell co-stimulation and Type II IFN Response (Figure 5E). Moreover, the expression of immune checkpoint molecules also differed significantly between the two risk groups, such as HHLA2, CD44 and TNFSF9 in the high-risk group were higher than those in the low-risk group, while other molecules were higher in the low-risk group (Figure 6A).

Finally, detailed Spearman correlation analysis was conducted. The risk score was positively correlated with common lymphoid progenitor, CD4⁺ T cells, Macrophage M1, Macrophage M0, Neutrophil, NK cells, Eosinophil and Monocyte. It was inversely associated with most other tumor-infiltrating immune cells (Figure 6B).

Wound healing assay, transwell assay and cell proliferation assay were performed in order to investigate the effect of LINP1 and SH3PXD2A-AS1 on invasion, migration, and proliferation of pancreatic cancer cells. As shown in Supplementary Figures S3A,B, LINP1 and SH3PXD2A-AS1 knockdown in PANC-1 inhibited the wound closure compared to the control cells. Similarly, LINP1 and SH3PXD2A-AS1 knockdown resulted in the decreased of invasive cells and cell viability (Supplementary Figures S3C–E) compared to the control group. These results indicated that downregulation of LINP1 and SH3PXD2A-AS1 can inhibit the proliferation, invasion and migration of pancreatic cancer cells.