Gene expression data and clinical information were downloaded from the TCGA-PAAD cohort FPKM file from the TCGA Information Portal (https://portal.gdc.cancer.gov/). After excluding 14 patients without RNA-seq or clinical information, 171 patients were included in this study. The 171 patients were divided into training cohort (n = 86) and testing cohort (n = 85). The training cohort were used for risk score model building, and the testing cohort and entire cohort were used for validation of the risk score model to test the robustness of the risk model.

TCGA data is open to the community and our research is based on the TCGA Data Access Policy and Publication Guidelines.

The ferroptosis-related mRNAs were obtained from the FerrDb database (http://www.zhounan.org/ferrdb), which comprised a total of 239 genes. Then, Pearson’s correlation analysis was performed for FRLs identified using the “limma” package. LncRNAs with R2 > .6 and p < .001 were considered to be FRLs.

First, the training cohort was used to construct a prognostic FRLs signature. Co-expressed FRLs were tested using the univariate Cox analysis to determine the prognostic value. Subsequently, based on the candidate lncRNAs with p < .01 in the univariate Cox screen analysis, a least absolute shrinkage and selection operator (LASSO) regression model was constructed by “glmnet” package. Morever, the resulting lncRNAs were introduced into a multivariate Cox model to obtain hazard ratios (HR) and regression coefficients for each lncRNA used in the construction of the final prognostic FRLs signature using “survival” package. Risk scores were calculated using the following equation. The risk score was calculated using the following formula: Risk   score = E x p lncRNA 1 × Coe f lncRNA 1 + E x p lncRNA 2 × Coe f lncRNA 2 + … + E x p lncRNAi × Coe f lncRNAi where E x p l n c R N A i is the expression of the ith selected lncRNA, and Coef l n c R N A i is its regression coefficient.

Subsequently, patients were divided into high- and low-risk groups based on the median risk score. The “survival” package was used to perform the Kaplan Meier (KM) survival and receiver operating characteristic (ROC) curve analysis of the performance of prognostic factors in terms of overall survival (OS). The testing cohort and the entire cohort were applied to validate the prognosis signature constructed based on the training cohort.

The KM analysis was performed in high-risk and low-risk groups to elucidate differences between subgroups based on gender (male and female), age (<65 and ≥65 years), grade (G1, G2 and G3, G4), T stage (T1, T2 and T3, T4) and N stage (N0 and N1). Univariate and multivariate Cox regression analysis of clinical characteristics and risk scores were performed to assess whether or not the risk score was an independent prognostic factor as well as to identify clinical characteristics that were independent prognostic factors.

In the classification of clinicopathological characteristics, the age, sex, grade, American Joint Committee on Cancer (AJCC) stage, T-stage, and N-stage were used to elucidate the independent risk factors associated with the prognosis. However, the M-stage was excluded from the study because of the high rate of missing data on the M-stage.

A nomogram was used to predict the 12-, 18-, and 24-month survival. The stability of the model was assessed using the consistency index (C-index) and calibration curves.

LncRNA-mRNA co-expression pairs with R ² > .4 and p < 0.05 were considered to be potential regulatory pathways. Subsequently, the co-expression network was visualized by Cytoscape software (version 3.8.2, http://www.cytoscape.org/).

To elucidate potential biological functions of FRLs, Gene Ontology (GO) (Gene Ontology Consortium, 2015) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) was performed on lncRNA-mRNAs. Then, significantly related biological functions and pathways, which may be potential pathways for ferroptosis-related lncRNA regulation, were listed.

Conversion of gene names and merging of the patient gene expression data with the clinical data were performed using PERL (version 5.30.2, http://www.perl.org/). All statistical analyses were performed in the R software (version 4.1.0).

Data from the TCGA-PAAD cohort were downloaded and processed, as reported previously. Inclusion criteria were: 1) a histological diagnosis of pancreatic adenocarcinoma; 2) available data on gene expression and clinicopathological characteristics; and 3) complete survival information and a follow-up duration >30 days. When a patient had two or more samples then the first sample was selected. After applying inclusion and exclusion criteria, 171 patients were enrolled. Table 1 summarizes the basic patient profile. Then, the 171 PAAD patients were randomly divided into a training cohort (n = 86) and a testing cohort (n = 85) at a 1:1 ratio. The flow chart of the study is shown in Figure 1.

First, 13,830 lncRNAs were identified by analyzing gene expression data from TCGA-PAAD. Altogether 239 ferroptosis-related genes (Supplementary Table S1) were expressed in patients with PAAD. Finally, 435 FRLs were selected. (R ² > 0.6 and p < .001).

The univariate Cox regression performed on 86 training PAAD samples with OS revealed 193 FRLs to be relevant with the PAAD prognosis (p < .01, Supplementary Table S2). Subsequently, LASSO regression of these genes, and eight FRLs were identified after LASSO regression (Figures 2A,B). Finally, multivariate Cox regression revealed six FRLs (LINC01705, AC068620.2, TRAF3IP2-AS1, AC092171.2, AC099850.3, and MIR193BHG, Figure 2C). A FRLs signature was created (Table 2) with the following risk score formula: Risk   score = ( ExpLINC 01705 × .0892 ) − ( ExpAC 068620.2 * 1.5521 ) − ( ExpTRAF 3 IP 2 − AS 1 * 2.7649 ) − ( ExpAC 092171.2 * 0.3138 ) + ( ExpAC 099850.3 * 0.1522 ) + ( ExpMIR 193 BHG ∗ 0.3023 )

Table 3 summarizes previous reports on the six FRLs. To the best of our knowledge, two lncRNAs (AC068620.2 and AC092171.2) have not been previously reported.

LINC01705, AC099850.3, and MIR193BHG were risk factors for PAAD, while three lncRNAs (AC068620.2, TRAF3IP2-AS1, and AC092171.2) were protective factors for PAAD (Figures 3A–F). The risk score formula applied to gene expression data of 86 training patients revealed a median risk score of 1.231. Accordingly, all cases could be classified as low-risk (n = 43) or high-risk (n = 43) (Figure 4A). The increased risk score was accompanied by gradually decreasing survival time and increased mortality (Figure 4B). Finally, the expression level of OS-related FRLs was presented in the form of heatmaps (Figure 4C). Regrettably, none of the patients in the high-risk group survived longer than 4 years, whereas many of those in the low-risk group survived longer than 7 years. After log-rank test, there was a significant difference in OS between the two groups (p = 1.135e-08, Figure 4D). Time-dependent ROC curves showed prediction of lncRNA biomarkers in the prognosis, with AUC of 0.786, respectively (Figure 4E). We then validated the a prognostic FRLs signature on the testing cohort and the entire cohort (Figures 4F–O). Survival analysis showed, significant differences in survival time between both the testing cohort and the entire cohort of patients, with significantly higher OS in the low-risk group than in the high-risk group. In addition, the AUCs for the testing cohort and the entire cohort were .723, .753, respectively, showing that the six FRLs prognostic signature has good robustness.

Univariate Cox, regression revealed the risk score to be a risk factor significantly associated with the prognosis of patients with PAAD (95% confidence interval (CI): 1.228–1.445, p < .001), with a higher risk score suggesting a worse prognosis. In addition, the N-stage (95% CI: 1.095–2.659, p = .018) and patient age (95% CI: 1.004–1.049, p = .022) were also closely associated with the prognosis. (Figure 5A). After controlling for clinical characteristics, risk score remained an independent indicator (hazard ratio = 1.372, 95% CI: 1.250–1.506, p < .001, Figure 5B, Table 4). However, the N-stage was not significantly associated with the prognosis in multivariate Cox regression. Moreover, higher accuracy of the risk score in predicting patients compared to other clinicopathological features was also confirmed in the ROC analysis. The prognostic value of FRLs (.771) was higher than that of age (.678), sex (.537), grade (.604), the AJCC stage (.581), the T-stage (.555), and the N-stage (.653; Figure 5C). Similarly, across all age, sex, T-stage, grade, and AJCC stage subgroups, patients with PAAD showed a lower OS rate in the high-risk group than in the low-risk group (Figure 6). All data suggested that the prognostic FRLs signature is significantly associated with the prognosis of patients with PAAD.

The nomogram, a common prognostic visualization tool in oncology, allows for the quantification of patient survival based on the inclusion of an index score. Based on multivariate Cox regression, the indicators, age and risk score, were included in the construction of the nomogram. Total scores were obtained by adding the individual indicator scores for age and risk score and predicting the probability of survival at 12, 18, and 24 months (Figure 7A). The nomogram had good stability at a consistency index of .697. Furthermore, calibration curves showed that predicted survival times at 12, 18, and 24 months were consistent compared to the reference line (Figures 7B–D). Furthermore, Table 5 lists other studies in which lncRNA predicted the prognosis of patients with PAAD.

The co-expression analysis of FRLs revealed a total of 61 mRNAs involved in formation of the co-expression network of 76 lncRNA-mRNA targets with R ² > .4 and p < .05 (Figure 8A; Supplementary Table S3). The GO analysis showed that it was involved in 476 biological processes, such as cellular responses to chemical stress and responses to oxidative stress, constituted 25 cellular components, including the phagosome assembly site and the apical portion of cells, and mediated 35 molecular functions, such as ubiquitin protein ligase binding and ubiquitin-like protein ligase binding (Figure 8B; Table 6). KEGG showed ferroptosis to be a significantly enriched pathway (Figure 8C; Table 7).