Integrated analysis of necroptosis-related lncRNAs for prognosis and immunotherapy of patients with pancreatic adenocarcinoma

Accumulating studies have revealed that necroptosis plays a vital role in the occurrence and development of pancreatic adenocarcinoma (PAAD). We aimed to construct a prognostic model for PAAD on the basis of necroptosis-related lncRNAs (NRLs). A coexpression network between necroptosis-related mRNAs and NRLs based on The Cancer Genome Atlas (TCGA) was constructed. Then, differentially expressed necroptosis-related lncRNAs (DENRLs) were screened from TCGA and Genotype-Tissue Expression project (GTEx) datasets. Univariate Cox regression (uni-Cox) analysis was performed on these DENRLs to identify lncRNAs significantly correlated with prognosis. Least absolute shrinkage and selection operator (LASSO) regression was performed for preventing overfitting on these lncRNAs. Multivariate Cox analysis (multi-Cox) was performed to establish a risk model based on lncRNAs that served as an independent prognostic factor. Next, the Kaplan–Meier analysis, time-dependent receiver operating characteristics (ROC), uni-Cox, multi-Cox regression, nomogram, and calibration curves were constructed to support the accuracy of the model. Gene set enrichment analysis (GSEA) and single-sample GSEA (ssGSEA) were also performed on risk groups, and it was found that the low-risk group was closely correlated with immune infiltration and immunotherapy. To further evaluate the immune differences between different clusters, we divided the patients into two clusters. Cluster 2 was more significantly infiltrated with immune cells and had higher immune scores. These results shed new light on the pathogenesis of PAAD based on NRLs and develop a prognostic model for diagnosing and guiding personalized immunotherapy of PAAD patients.

Introduction

Pancreatic adenocarcinoma (PAAD) is the fourth leading cause of cancer-related deaths worldwide, with 5-year survival rates of less than 11% (Siegel et al., 2022). Although significant progress has been made in the treatment of pancreatic cancer in the past few decades, the postoperative prognosis of pancreatic cancer patients is not ideal because of the insidious onset and high recurrence and metastasis rate of pancreatic cancer (Neoptolemos et al., 2018; Mizrahi et al., 2020). Compared with a single clinical marker, a model integrating multiple markers can more accurately predict patient prognosis and provide targets for therapy on the basis of the development of large-scale genome sequencing technologies. Therefore, it is important to construct a molecular model that can predict patient prognosis and guide individualized therapy.

Necroptosis is a form of programmed cell death that can serve as a mode for cells to overcome apoptotic resistance (Gong et al., 2019), which plays an important role in inflammation, microbial infection, and tumor metastasis (Seifert et al., 2016; Place et al., 2021; Liu Z. et al., 2021; Yan et al., 2022). Long noncoding RNAs (lncRNAs), which are defined as RNA molecules more than 200 nucleotides in length and have limited protein-coding ability, play an important part in the development and progression of different cancers, including PAAD (Liu et al., 2019; Zhou et al., 2020). Accumulating evidence revealed that lncRNAs are involved in the occurrence and development of necroptosis (Wang et al., 2020; Harari-Steinfeld et al., 2021). However, the role of prognosis and potential molecular mechanism associated with NRLs is currently unclear in PAAD.

In the present study, we established a prognostic model based on NRLs and found that it was able to divide patients into two groups with significant differences in prognosis and immune infiltration, which may provide diverse insights into a clinical therapy for PAAD patients.

Materials and methods

Acquisition of datasets and patients’ information of PAAD

The PAAD RNA-seq dataset and GTEx data were downloaded from UCSC Xena (http://xena.ucsc.edu/) (including 178 tumor samples and 171 normal samples). Corresponding clinical data (including age, gender, survival time, survival status, grade, and pathological stage) were also obtained from UCSC Xena. Then we converted the data of GTEx in the form of log2 (X+0.001) to log2 (X+1) to be consistent with the data of PAAD by Strawberry Perl. The study excluded patients with missing survival information to reduce statistical bias in this analysis.

Identification of DENRLs

A total of 67 necroptosis-related genes (NRGs) were collected from previous literature (Zhao et al., 2021). Then, 356 necroptosis-related lncRNAs (NRLs) were obtained from a correlation analysis performed between NRGs and lncRNAs of PAAD (|Pearson correlation coefficients| > 0.7, and p < 0.001). The mRNA–lncRNA coexpression network was visualized by the “igraph” R package. The Wilcox test was performed to explore the expression of 356 NRLs between PAAD and normal tissues with the “limma” package (Ritchie et al., 2015) in R (Rx64 4.1.2, https://www.rstudio.com/) (Log2 fold change (FC) > 1, false discovery rate (FDR) < 0.05). Finally, 172 DENRLs were obtained and visualized by the “pheatmap” R package.

Construction and validation of a prognostic NRL model for PAAD

DENRLs’ expression data merged with the survival information of 176 patients were retrieved and regrouped into the train, test, and entire group randomly by the Perl script and the “caret” R package with a ratio of 1:1. Based on uni-Cox analysis with p < 0.01 as the threshold, prognosis-related lncRNAs were screened and incorporated into a Lasso regression analysis to conduct multivariate Cox proportional hazard regression analysis with p < 0.05 as the threshold. After narrowing the range of prognosis-related lncRNAs and preventing overfitting, a model was constructed by using the R package “glmnet” (Friedman et al., 2010; Simon et al., 2011). The risk score is calculated by the following formula:

$\text{risk}\text{score}={\displaystyle \sum }_{i=1}^n\text{coefficient}\left({\text{lncRNA}}^i\right)\ast \text{expression}\left({\text{lncRNA}}^i\right).$

The patients were divided into low-risk and high-risk groups by the median risk score. The Chi-square test was used to evaluate the correlation of clinical factors with the model and visualized by the R package “survminer.” Univariate Cox and multivariate Cox regression analyses were developed to determine whether the risk score and clinical characteristics were independent predictors. The R package “survivalROC” was used to evaluate the accuracy of the constructed model. The larger the AUC (area under the ROC curve), the higher the prediction accuracy of the model.

Construction and validation of the nomogram

To improve the predictive accuracy of NRLs in the model of PAAD patients, a nomogram was established for forecasting 1-, 2-, and 3-year OS integrating risk scores, age, gender, grade, T stage, N stage, and pathological stage by the “rms” R package. Then, the calibration curve was used to illustrate the discrimination between the nomogram-predicted value and the true value.

Gene set enrichment analysis between high- and low-risk groups

We used gene set enrichment analysis (GSEA) software 4.2.1 (http://www.gsea-msigdb.org/gsea/index.jsp) to investigate the significantly predefined biological pathways between the high- and low-risk groups with the threshold: p < 0.01 and FDR <0.1. c2. cp.kegg.v7.1. symbols.gmt was chosen from the Molecular Signatures Database (MSigDB http://software.broadinstitute.org/gsea/msigdb/index.jsp) as the reference file. All these results were visualized by ggplot2, grid, and gridExtra R packages.

Exploration of immune infiltration and immune checkpoints

To figure out the correlation between NRL signatures in the model and the tumor immune microenvironment, we calculated the infiltration values for PAAD patients based on seven algorithms: XCELL (Aran et al., 2017), TIMER (Li et al., 2017), QUANTISEQ (Finotello et al., 2019), MCPCOUNTER (Becht et al., 2016), EPIC (Racle et al., 2017), CIBERSORT−ABS (Newman et al., 2015), and CIBERSORT (Newman et al., 2015). Then, Spearman correlation analysis visualized in a bubble chart was used to evaluate the relationship between the immune and risk scores. Subsequently, single-sample GSEA (ssGSEA) was performed for investigating the difference between immune cells and immune function between high- and low-risk groups. The stromal score, immune score, and ESTIMATE score of each patient in different risk groups were calculated. The differences were visualized by the “ggpubr” R package. In addition, different expressions of immune checkpoint genes between the risk groups were evaluated by the “ggpubr” R package.

Clusters based on the expression of lncRNAs in the model

As the literature illustrated previously, different subtypes usually contributed to different immune microenvironments (DeBerardinis, 2020). To explore the immunotherapy response of different subtypes, subgroups based on the expression of lncRNAs in the model by ConsensusClusterPlus (CC) (Wilkerson and Hayes, 2010) R package were performed to investigate the response to immunotherapy within them. Parameters of the consensus clustering model were set as maxK = 9, reps = 50, pItem = 0.8, pFeature = 1, clusterAlg = "km,” and distance = "euclidean.” Principal component analysis (PCA), T-distributed stochastic neighbor embedding (t-SNE), and Kaplan–Meier survival curve were constructed by the “Rtsne” R package.

Cell culture and quantitative real-time RT-PCR (RT-qPCR)

Normal human pancreatic ductal epithelial cells hTERT-HPNE and two pancreatic cancer cells MiaPaCa-2 and BxPC-3 were purchased from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). HPNE and MiaPaCa-2 were cultured in Dulbecco’s Modified Eagle Medium (Gibco, United States) with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. BxPC-3 was cultured in Roswell Park Memorial Institute (RPMI) medium (Gibco, United States) with 10% FBS. A humidified atmosphere containing 5% CO₂ at 37 °C was provided to culture all cells. A Fastgen 200 RNA isolation kit (Shanghai, China) was used to extract the total RNA according to the protocol provided. Then, the total RNA was reverse transcribed into cDNA by using the Prime Script RT reagent kit (TaKaRa, Dalian, China). RT-qPCR was performed using the CFX Manager 2.1 fluorescent quantitative PCR kit (Bio-Rad Laboratories, United States). The primer sequences were as follows: linc01089 (human), sense 5′-CCT​CTA​GCA​GAG​TGC​CTT​GG-3′ and antisense 5′-AGG​TAA​CCG​GGG​TCA​GAT​CA-3’; GAPDH (human), sense 5′-GGA​GCG​AGA​TCC​CTC​CAA​AAT-3′ and antisense 5′-GGC​TGT​TGT​CAT​ACT​TCT​CAT​GG-3’. GAPDH was used as an internal control.

Results

Identification of necroptosis-related lncRNAs

A total of 178 tumor samples and 171 normal samples (167 samples from GTEx) were obtained from UCSC Xena based on TCGA and GTEx. Then we identified 356 NRLs and expression data by Pearson correlation analysis based on 69 NRGs (|Pearson correlation coefficients| > 0.7, and p < 0.001). The network between NRGs and NRLs are shown in Figure 1A. Subsequently, 172 DENRLs were screened between tumor and normal samples, of which 66 were downregulated and 106 were upregulated (Log2 fold change (FC) > 1, false discovery rate (FDR) < 0.05, Figure 1B).

FIGURE 1

FIGURE 1. Identification and selection of NRL prognostic signatures in PAAD. (A) Coexpression network between NRGs and NRLs (|Pearson correlation coefficients| > 0.7, and p < 0.001). (B) Volcano plot of 172 DENRLs. (C) The expression profiles of 36 prognostic-related NRLs. (D) The Sankey diagram of prognostic NRLs, NRGs, and risk type. (E) The LASSO coefficient profile of necroptosis-related lncRNAs. (F) The 10-fold cross-validation for variable selection in the LASSO model.

Construction and validation of the prognostic model based on TCGA cohorts

To predict the prognosis and provide therapeutic targets for PAAD patients, we constructed a prognostic model based on NRLs. A total of 36 prognostic-related NRLs were screened by univariate Cox regression (p < 0.01) and visualized with a heatmap based on their expression between tumor and normal samples (Figure 1C). The correlations between prognostic NRLs, NRGs, and the risk type were displayed by a Sankey diagram (Figure 1D). To narrow the range of prognosis-related lncRNAs and avoid overfitting, we performed Lasso regression and extracted nine NRLs (Figures 1E,F). Then, multivariate Cox proportional hazard regression analysis was conducted and five NRLs were screened to establish a prognostic NRL signature model (U62317.1, AC099850.3, LINC01089, LINC01133, and LINC00996). The risk score formula was as follows: (0.4675 × U62317.1 expression) + (0.3947 × AC099850.3 expression) + (−0.5810 × LINC01089 expression) + (0.2759 × LINC01133 expression) + (−2.0985 × LINC00996 expression). With this formula, we calculated the risk score for each patient based on personalized NRL expression levels in train, test, and entire groups and grouped them into low- and high-risk groups by the median risk score. All these results implied a significant statistical difference in the OS between the high- and low-risk groups, indicating that these developed signatures effectively predict prognosis (Figures 2A–C, Supplementary Figure S1). Subsequently, univariate and multivariate Cox regression analyses were performed to determine whether risk scores and clinical characteristics are independent prognostic factors for the OS in patients with PAAD. The hazard ratio (HR) (95% CI) for the risk score was 1.272 (1.189–1.361) in univariate Cox regression analysis (p < 0.001, Figure 2D) and 1.262 (1.175–1.355) in multivariate Cox regression analysis (p < 0.001, Figure 2E), which indicates that the NRLs in our model are independent prognostic indicators. In addition, the model was evaluated by ROC analysis at 1, 3, and 5 years to verify its predictive value, and the corresponding AUC was 0.759, 0.766, and 0.783, respectively. (Figure 2F). Furthermore, we combined the clinical information of age, gender, grade, and stage to construct a time-dependent ROC curve to predict the 1-year survival rate of patients with PAAD, and the AUC was 0.540, 0.565, 0.597, and 0.483 (Figure 2G). To evaluate the importance of NRLs in the progression of PAAD, we investigated the correlation between the risk score and clinical characteristics. As shown in Figure 3, the OS rate in the high-risk group of subgroups separated by age, gender, T stage, N stage, M stage, clinical stage, and grade was lower than that of the low-risk group (p < 0.05). Taken together, the risk model we constructed was feasible and reliable in identifying the risk of PAAD patients.

FIGURE 2

FIGURE 2. Prognosis value of the five NRL models in the train, test, and entire sets. (A) Risk score distribution, scatter plots, and KM survival curve analysis of patients with PAAD in the training cohort. (B) Risk score distribution, scatter plots, and KM survival curve analysis of patients with PAAD in the testing cohort. (C) Risk score distribution, scatter plots, and KM survival curve analysis of patients with PAAD in the entire cohort. (D) Univariate Cox regression showing that the age, grade, and risk score were associated with the OS (p < 0.05). (E) Multivariate Cox regression showing that the risk score (p < 0.001) was an independent prognostic indicator of the OS in patients with PAAD. (F) 1-, 3-, and 5-year ROC curves of the entire set. (G) 1-year ROC curves of the risk score and clinical characteristics.

FIGURE 3

FIGURE 3. Survival rate of PAAD patients with high- and low-risk patients in the subgroups based on clinicopathological characteristics.

Construction and verification of the nomogram

To further enhance the predictive power of the risk model, we constructed a nomogram integrated with clinicopathological factors to predict 1-, 2-, and 3-year OS (Figure 4A). Based on the nomogram score, the 1-, 2-, and 3-year OS could be well predicted. In addition, as shown in Figure 4B, good consistency between the nomogram-predicted OS and observed OS could be found at 1, 2, and 3 years with calibration curves. The aforementioned results indicated that the nomogram had a good predictive value in PAAD patients.

FIGURE 4

FIGURE 4. Construction and validation of the nomogram. (A) Prognostic nomogram that integrated risk score, age, gender, grade, T, N, and tumor stage predicted the probability of 1-, 2-, and 3-year OS. (B) Calibration curves for 1-, 2-, and 3-year OS.

Gene set enrichment analysis between the high- and low-risk groups

To further investigate the biological functions of the prognostic signature, we performed GSEA between the risk groups in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway in the entire set (Figure 5A). The P53 signaling pathway (NES = 1.95, p = 0.000), cell cycle (NES = 1.91, p = 0.006), base excision repair (NES = 1.90, p = 0.000), pyrimidine metabolism (NES = 1.86, p = 0.000), and homologous recombination (NES = 1.85, p = 0.002) were enriched in the high-risk group. Neuroactive ligand-receptor interactions (NES = -1.78, p = 0.002), adipocytokine signaling pathway (NES = -1.75, p = 0.004), FC epsilon RI signaling pathway (NES = -1.75, p = 0.008), calcium signaling pathway (NES = -1.74, p = 0.008), and type II diabetes mellitus (NES = -1.70, p = 0.006) were enriched in PAAD samples with the low-risk group. These results indicated that the prognostic signature in the model may be associated with tumor progression and perineural invasion.

FIGURE 5

FIGURE 5. Investigation of tumor immune factors. (A) Top five pathways significantly enriched in the high- and low-risk groups were exhibited by GSEA. (B) Immune cell bubble of risk groups. (C) Comparison of immune-related scores between low- and high-risk groups. (D) Difference of immune functions in risk groups. (E) Difference of immune cells in risk groups. (F) Differential expression of 29 checkpoints’ expression in risk groups.

Investigation of the correlation of immune infiltration in risk groups

A stronger correlation existed between the various immune cells and the low-risk group on different platforms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT) exhibited with an immune cell bubble chart (Figure 5B). Then we found that the low-risk group had a higher immune score and ESTIMATE score, displaying a significant difference from the high-risk group (Figure 5C). In addition, the stronger correlation between immune function-associated biomarkers and ssGSEA scores showed that the low-risk group had a higher immune infiltration status (Figures 5D,E). The expression levels of immune checkpoints and/or their ligands may reflect the effect of immunotherapy. Then the correlation between the expression of 29 immune checkpoint genes and risk scores was investigated (Figure 5F). Low-risk patients exhibited a higher expression of 25 immune checkpoint genes, namely, BTLA, CD200, NRP1, LAIR1, CD244, LAG3, CD200R1, ICOS, CD40LG, CTLA4, CD48, CD28, HAVCR2, ADORA2A, KIR3DL1, TNFSF14, CD160, IDO2, IDO1, PDCD1, TNFRSF8, CD27, BTNL2, TIGIT, and CD86, while the remaining four were highly expressed in the high-risk group. Collectively, these results indicated that NRL signatures in the model affect tumor-infiltrating immune cell (TIC) infiltration and can be a candidate biomarker for immunotherapy.

Clusters

To evaluate the immune differences between different clusters, patients were divided into different clusters by the “ConsensusClusterPlus” R package based on the expression of NRLs in the model (Supplementary Figure S2). According to the results of clustering analysis, the optimal number of clusters was 2 (Figure 6A). The correlation between two clusters and risk groups is vividly shown in Figure 6B by a Sankey diagram. T-SNE and PCA indicated that the risk groups and two clusters could be distinguished clearly (Figure 6D). Furthermore, cluster 2 (C2) had better OS than cluster 1 (C1) (Figure 6C). As shown in Figures 6F,G, C2 had a higher immune score, estimate score, and immune cell infiltration, signifying a higher immune infiltration status than C1. With the immune checkpoint gene comparison, we found that C2 tends to have a higher expression of 19 immune checkpoint genes, namely, BTLA, CD200, NRP1, LAIR1, CD200R1, ICOS, CD40LG, CTLA4, CD48, CD28, ADORA2A, KIR3DL1, TNFSF14, CD160, IDO2, PDCD1, CD27, TIGIT, and CD86, while the remaining six immune checkpoint genes (TNFRSF14, LGALS9, TNFSF9, TNFSF15, HHLA2, and CD70) were highly expressed in C1 (Figure 6E). Based on the NRLs of different clusters, we could further explore the correlation between immunotherapy and NRL signatures in PAAD patients.

FIGURE 6

FIGURE 6. Distinction between different clusters. (A) Patients divided into two clusters by ConsensusClusterPlus. (B) Sankey diagram of different clusters and risk groups. (C) Kaplan–Meier survival curves of OS in clusters. (D) PCA and t-SNE of risk groups and clusters. (E) Differential expression of 25 checkpoint genes in different clusters. (F) Comparison of immune-related scores between clusters 1 and 2. (G) Heatmap of immune cells in clusters.

Validation of the expression and prognostic value of linc01089 in PAAD

The expression and prognosis value of linc01089 was explored based on RT-qPCR and TCGA dataset. Consistent with the previous results (Figure 1C), the expression of linc01089 was lower in PAAD cells than in HPNE (Figure 7A). Furthermore, the expression of linc01089 in stages III and IIV was lower than that in stages I and II (Figure 7B). Univariate and multivariate analyses revealed that the N stage and linc01089 expression were independent factors affecting the prognosis of PAAD patients (Figures 7C,D).

FIGURE 7

FIGURE 7. Validation of the expression and clinical characteristics of linc01089 in PAAD. (A) Expression of linc01089 was downregulated in PAAD cells than that in HPNE. (B) Correlation between linc01089 and the pathology stage in PAAD. (C,D) Univariate and multivariate analyses indicated that linc01089 and N stage were independent prognostic factors affecting the OS of PAAD patients.

Discussion

Pancreatic cancer is a highly malignant tumor of the digestive system. Because of its insidious symptoms and rapid disease progression, most patients present with locally advanced or metastatic disease at diagnosis and the treatment effect is not ideal (Park et al., 2021). Currently, the gene signatures of distinct molecular subtypes based on cell activities including autophagy (Hu et al., 2020), metabolic (Yang et al., 2020), immune (Kandimalla et al., 2020), and so on, show good predictive value for individual risk assessment in cancers. Given this reason, constructing a molecular prediction model in PAAD for guiding personalized therapy and predicting the prognosis is particularly important.

Necroptosis is a form of cell death which occurs downstream of protein kinase C-related kinase 1 and 3 (PRK1 and RIPK3) as members of the oligomeric complex termed the necrosome (Galluzzi et al., 2012). Accumulating studies illustrated that necroptosis plays a tumor inhibitory role (Newton, 2015) and reduces the incidence of metastasis by the accumulation of high ROS levels (Lawlor et al., 2015). LncRNAs involved in necroptosis have been continuously identified (Jiang et al., 2021). Harari-Steinfeld et al. (2021) found that lncRNA H19-derived microRNA-675 promotes liver necroptosis by targeting FADD. Wang et al. (2020) found that ammonia regulates chicken tracheal cell necroptosis via the lncRNA-107053293/miR-148a-3p/FAF1 axis. Although several NRL signatures of cancer have been reported (Wang and Liu, 2021; Zhao et al., 2021), there is no study associated with NRL signatures of pancreatic cancer. In this study, we established a prognostic model with five NRLs (U62317.1, AC099850.3, LINC01089, LINC01133, and LINC00996). On the basis of univariate and multivariate Cox regression analyses, the NRL signatures in the model are independent prognostic indicators, implying a moderate predictive value. In addition, the ROC curve and nomogram displayed a good predictive effect to predict prognosis and guide immunotherapy.

Among these NRLs in the model, all were associated with tumor progression in previous research. U62317.1 is identified as an lncRNA biomarker for the early diagnosis of oral squamous cell carcinoma (Li Y. et al., 2020). AC099850.3 is a member of the biomarkers that participate in autophagy, immune, and stemness (Jia et al., 2020; Wu et al., 2021; Zhang et al., 2021). LINC01089 inhibits the progression of gastric cancer (Guo and Li, 2020), cervical cancer (Li et al., 2021), colorectal cancer (Li and Guo, 2020), and lung cancer (Li X. et al., 2020) by sponging miRNA. LINC01133 aggravates the epithelial-mesenchymal transition through the Wnt/β-catenin pathway by silencing AXIN2 (Liu Y. et al., 2021). LINC00996 is a potential therapeutic target in pulmonary adenocarcinoma and squamous cell carcinoma (Yan et al., 2021). However, studies on cell activities and NRLs are still lacking. Thus, more experimental evidence of correlations between NRLs and PAAD is needed.

Growing evidence suggests that immune cells contribute to tumor progression when present in the tumor microenvironment (TME) (Hinshaw and Shevde, 2019). The efficacy of immunotherapy and overall survival were significantly correlated with the composition of the TME (Petitprez et al., 2020). In this study, patients with a low-risk score had a higher infiltration level with immune cells including B cells, CD8⁺T cells, dendritic cells, and NK cells and a stronger correlation with immune function, confirming the function of NRLs in tumor immune infiltration. Immune checkpoint blockade (ICB) enhances the efficiency of anti-tumoral immune response by the way of immune infiltration (Petitprez et al., 2020). Interestingly, CTLA4, IDO1, and PDCD1 were highly expressed in the low-risk group, implying that low-risk patients may benefit more from immunotherapy. There is a difference between different subtypes and immune scores which contribute to different prognoses and immunotherapy responses (DeBerardinis, 2020). Then we divided patients into two clusters on the basis of NRL expressions in the model. As expected, differences existed between clusters in the immune score and immune cell infiltration. C2 had a higher immune score and estimate score. At the same time, C2 had more immune cell infiltration and a higher expression of CTLA4 and PDCD1, which confirmed that C2 may benefit more from immunotherapy than C1. Above all, these NRLs could accurately predict the prognosis of PAAD patients and guide personalized immunotherapy.

There are still several limitations in the current study. First, external validation with data integrated with complete clinical and lncRNA information will be searched and analyzed subsequently. Second, functional experiments including cell and animal experiments will be in our plans to explore molecular mechanisms in-depth.

Conclusion

All in all, we constructed a model with NRLs to predict prognosis and guide personalized immunotherapy in PAAD patients. Targeting these lncRNAs will be a promising way for systemic therapy failure and new pathways of immunotherapy. Therefore, molecular mechanisms between NRLs, necroptosis, and PAAD are worth investigating.

Data availability statement

PAAD RNA-seq dataset and Genotype-Tissue Expression project (GTEx) data used in this study are available via the UCSC Xena (http://xena.ucsc.edu/); further inquiries can be directed to the corresponding author.

Author contributions

JM and ZC contributed equally to this work and share first authorship. QM, ZW, and ZW contributed to the conception and design of the study. QW and WZ performed the statistical analysis. JM and ZC wrote the first draft of the manuscript. JL, SW, WQ, and CZ performed reviewing and editing. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This study was funded by Innovation Capability Support Program of Shaanxi (2022 PT-35), Institutional Foundation of The First Affiliated Hospital of Xi’an Jiaotong University (2021QN-24), and Natural Science Basic Research Program of Shaanxi (2020JQ-510).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.940794/full#supplementary-material

References

Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18 (1), 220. doi:10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

DeBerardinis, R. J. (2020). Tumor microenvironment, metabolism, and immunotherapy. N. Engl. J. Med. 382 (9), 869–871. doi:10.1056/NEJMcibr1914890

PubMed Abstract | CrossRef Full Text | Google Scholar

Finotello, F., Mayer, C., Plattner, C., Laschober, G., Rieder, D., Hackl, H., et al. (2019). Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 11 (1), 34. doi:10.1186/s13073-019-0638-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

Galluzzi, L., Vitale, I., Abrams, J. M., Alnemri, E. S., Baehrecke, E. H., Blagosklonny, M. V., et al. (2012). Molecular definitions of cell death subroutines: recommendations of the nomenclature committee on cell death 2012. Cell Death Differ. 19 (1), 107–120. doi:10.1038/cdd.2011.96

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, Y., Fan, Z., Luo, G., Yang, C., Huang, Q., Fan, K., et al. (2019). The role of necroptosis in cancer biology and therapy. Mol. Cancer 18 (1), 100. doi:10.1186/s12943-019-1029-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., and Li, M. (2020). LINC01089 is a tumor-suppressive lncRNA in gastric cancer and it regulates miR-27a-3p/TET1 axis. Cancer Cell Int. 20, 507. doi:10.1186/s12935-020-01561-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Harari-Steinfeld, R., Gefen, M., Simerzin, A., Zorde-Khvalevsky, E., Rivkin, M., Ella, E., et al. (2021). The lncRNA H19-derived MicroRNA-675 promotes liver necroptosis by targeting FADD. Cancers (Basel) 13 (3), 411. doi:10.3390/cancers13030411

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinshaw, D. C., and Shevde, L. A. (2019). The tumor microenvironment innately modulates cancer progression. Cancer Res. 79 (18), 4557–4566. doi:10.1158/0008-5472.CAN-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, D., Jiang, L., Luo, S., Zhao, X., Hu, H., Zhao, G., et al. (2020). Development of an autophagy-related gene expression signature for prognosis prediction in prostate cancer patients. J. Transl. Med. 18 (1), 160. doi:10.1186/s12967-020-02323-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, Y., Chen, Y., and Liu, J. (2020). Prognosis-predictive signature and nomogram based on autophagy-related long non-coding RNAs for hepatocellular carcinoma. Front. Genet. 11, 608668. doi:10.3389/fgene.2020.608668

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, N., Zhang, X., Gu, X., Li, X., and Shang, L. (2021). Progress in understanding the role of lncRNA in programmed cell death. Cell Death Discov. 7 (1), 30. doi:10.1038/s41420-021-00407-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandimalla, R., Tomihara, H., Banwait, J. K., Yamamura, K., Singh, G., Baba, H., et al. (2020). A 15-gene immune, stromal, and proliferation gene signature that significantly associates with poor survival in patients with pancreatic ductal adenocarcinoma. Clin. Cancer Res. 26 (14), 3641–3648. doi:10.1158/1078-0432.CCR-19-4044

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawlor, K. E., Khan, N., Mildenhall, A., Gerlic, M., Croker, B. A., D'cruz, A. A., et al. (2015). RIPK3 promotes cell death and NLRP3 inflammasome activation in the absence of MLKL. Nat. Commun. 6, 6282. doi:10.1038/ncomms7282

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., and Guo, X. (2020). LINC01089 blocks the proliferation and metastasis of colorectal cancer cells via regulating miR-27b-3p/HOXA10 axis. Onco. Targets. Ther. 13, 8251–8260. doi:10.2147/OTT.S256148

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). Timer: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77 (21), e108–e110. doi:10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Lv, F., Li, F., Du, M., Liang, Y., Ju, S., et al. (2020a). LINC01089 inhibits tumorigenesis and epithelial-mesenchymal transition of non-small cell lung cancer via the miR-27a/SFRP1/Wnt/β-catenin Axis. Front. Oncol. 10, 532581. doi:10.3389/fonc.2020.532581

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Cao, X., and Li, H. (2020b). Identification and validation of novel long non-coding RNA biomarkers for early diagnosis of oral squamous cell carcinoma. Front. Bioeng. Biotechnol. 8, 256. doi:10.3389/fbioe.2020.00256

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Han, Y., Liang, X., and Zhao, M. (2021). LINC01089 inhibits the progression of cervical cancer via inhibiting miR-27a-3p and increasing BTG2. J. Gene Med. 23 (1), e3280. doi:10.1002/jgm.3280

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, M., Zhong, J., Zeng, Z., Huang, K., Ye, Z., Deng, S., et al. (2019). Hypoxia-induced feedback of HIF-1α and lncRNA-CF129 contributes to pancreatic cancer progression through stabilization of p53 protein. Theranostics 9 (16), 4795–4810. doi:10.7150/thno.30988

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Tang, T., Yang, X., Qin, P., Wang, P., Zhang, H., et al. (2021a). Tumor-derived exosomal long noncoding RNA LINC01133, regulated by Periostin, contributes to pancreatic ductal adenocarcinoma epithelial-mesenchymal transition through the Wnt/β-catenin pathway by silencing AXIN2. Oncogene 40 (17), 3164–3179. doi:10.1038/s41388-021-01762-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Nailwal, H., Rector, J., Rahman, M. M., Sam, R., Mcfadden, G., et al. (2021b). A class of viral inducer of degradation of the necroptosis adaptor RIPK3 regulates virus-induced inflammation. Immunity 54 (2), 247–258. e247. doi:10.1016/j.immuni.2020.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizrahi, J. D., Surana, R., Valle, J. W., and Shroff, R. T. (2020). Pancreatic cancer. Lancet 395 (10242), 2008–2020. doi:10.1016/S0140-6736(20)30974-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Neoptolemos, J. P., Kleeff, J., Michl, P., Costello, E., Greenhalf, W., Palmer, D. H., et al. (2018). Therapeutic developments in pancreatic cancer: current and future perspectives. Nat. Rev. Gastroenterol. Hepatol. 15 (6), 333–348. doi:10.1038/s41575-018-0005-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Newton, K. (2015). RIPK1 and RIPK3: critical regulators of inflammation and cell death. Trends Cell Biol. 25 (6), 347–353. doi:10.1016/j.tcb.2015.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, W., Chawla, A., and O'reilly, E. M. (2021). Pancreatic cancer: a review. JAMA 326 (9), 851–862. doi:10.1001/jama.2021.13027

PubMed Abstract | CrossRef Full Text | Google Scholar

Petitprez, F., Meylan, M., De Reynies, A., Sautes-Fridman, C., and Fridman, W. H. (2020). The tumor microenvironment in the response to immune checkpoint blockade therapies. Front. Immunol. 11, 784. doi:10.3389/fimmu.2020.00784

PubMed Abstract | CrossRef Full Text | Google Scholar

Place, D. E., Lee, S., and Kanneganti, T. D. (2021). PANoptosis in microbial infection. Curr. Opin. Microbiol. 59, 42–49. doi:10.1016/j.mib.2020.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Racle, J., De Jonge, K., Baumgaertner, P., Speiser, D. E., and Gfeller, D. (2017). Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife 6, e26476. doi:10.7554/eLife.26476

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Seifert, L., Werba, G., Tiwari, S., Giao Ly, N. N., Alothman, S., Alqunaibit, D., et al. (2016). The necrosome promotes pancreatic oncogenesis via CXCL1 and Mincle-induced immune suppression. Nature 532 (7598), 245–249. doi:10.1038/nature17403

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. CA. Cancer J. Clin. 72 (1), 7–33. doi:10.3322/caac.21708

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization paths for cox's proportional hazards model via coordinate descent. J. Stat. Softw. 39 (5), 1–13. doi:10.18637/jss.v039.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, N., and Liu, D. (2021). Identification and validation a necroptosisrelated prognostic signature and associated regulatory axis in stomach adenocarcinoma. Onco. Targets. Ther. 14, 5373–5383. doi:10.2147/OTT.S342613

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Shi, Q., Wang, S., Zhang, H., and Xu, S. (2020). Ammonia regulates chicken tracheal cell necroptosis via the LncRNA-107053293/MiR-148a-3p/FAF1 axis. J. Hazard. Mat. 386, 121626. doi:10.1016/j.jhazmat.2019.121626

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Wei, H., Liu, G., and Zhang, Y. (2021). Bioinformatics profiling of five immune-related lncRNAs for a prognostic model of hepatocellular carcinoma. Front. Oncol. 11, 667904. doi:10.3389/fonc.2021.667904

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, T., Ma, G., Wang, K., Liu, W., Zhong, W., Du, J., et al. (2021). The immune heterogeneity between pulmonary adenocarcinoma and squamous cell carcinoma: a comprehensive analysis based on lncRNA model. Front. Immunol. 12, 547333. doi:10.3389/fimmu.2021.547333

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, J., Wan, P., Choksi, S., and Liu, Z. G. (2022). Necroptosis and tumor progression. Trends Cancer 8 (1), 21–27. doi:10.1016/j.trecan.2021.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Huang, X., Liu, Z., Qin, W., and Wang, C. (2020). Metabolism-associated molecular classification of hepatocellular carcinoma. Mol. Oncol. 14 (4), 896–913. doi:10.1002/1878-0261.12639

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Cheng, M., Fan, Z., Jin, Q., Cao, P., Zhou, G., et al. (2021). Identification of cancer cell stemness-associated long noncoding RNAs for predicting prognosis of patients with hepatocellular carcinoma. DNA Cell Biol. 40 (8), 1087–1100. doi:10.1089/dna.2021.0282

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., Liu, H., Zhou, X., Fang, D., Ou, X., Ye, J., et al. (2021). Necroptosis-related lncRNAs: Predicting prognosis and the distinction between the cold and hot tumors in gastric cancer. J. Oncol. 2021, 6718443. doi:10.1155/2021/6718443

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Yi, C., Yi, Y., Qin, W., Yan, Y., Dong, X., et al. (2020). LncRNA PVT1 promotes gemcitabine resistance of pancreatic cancer via activating Wnt/β-catenin and autophagy pathway through modulating the miR-619-5p/Pygo2 and miR-619-5p/ATG14 axes. Mol. Cancer 19 (1), 118. doi:10.1186/s12943-020-01237-y

PubMed Abstract | CrossRef Full Text | Google Scholar