We downloaded the transcriptome RNA sequencing data of LUAD samples from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Meanwhile, we obtained the corresponding clinical parameters of these patients and excluded patients with missing overall survival (OS) or poor OS (less than 60 days) to reduce statistical bias in this analysis.

A list of 67 necroptosis genes was obtained from previously reported literature (Zhao et al., 2021). The correlations between 67 necroptosis-related genes and lncRNA expression were analyzed via Pearson correlation analysis. All NRlncRNAs (2,154) should conform to the standard of correlation coefficients (|Pearson R|) >0.4 and p <0.001. Then, we obtained 1,061 differentially expressed lncRNAs [log2 fold change > 1, false discovery rate (FDR) <0.05] after screening the synthetic data matrix by Strawberry Perl V-5.30.0 (https://www.perl.org/) and R software V-4.1.2 (https://www.r-project.org/) with limma R package.

The entire 481 TCGA set of LUAD was divided into a train risk set and a test risk set randomly by the caret R package. The ratio was 1:1. The train set was used to construct a necroptosis-related lncRNA signature, and the test set and entire set were applied to validate the signature. Combined with the clinical information of LUAD in TCGA, we screened and obtained 40 NRlncRNAs linked to OS significantly by univariate Cox (uni-Cox) regression analysis (p < 0.05). Subsequently, we performed least absolute shrinkage and selection operator (LASSO) Cox analysis (using the penalty parameter estimated by 10-fold cross-validation) via the glmnet R package to screen out optimal lncRNAs associated with LUAD prognosis. This method aims to prevent over-fitting during modeling. Finally, a prognostic risk signature based on the optimal lncRNAs was established with the multivariate Cox (multi-Cox) regression analysis, and the risk score of every patient with LUAD was calculated based on the following formula: risk score = ∑ i = 1 n Coef ( i ) × Expr ( i )

Coef(i) and Expr(i) represent the regression coefficient of the multi-Cox regression analysis for each lncRNA and each lncRNA expression level, respectively. The patients were stratified into low- and high-risk groups, with the risk score as the cutoff. Kaplan–Meier method and log-rank test were conducted to analyze whether there is a difference in the OS of LUAD patients between the low- and high-risk groups using the survival R package.

We evaluated the prognostic value of the established risk signature between the model and the clinical characteristics via chi-square test. Uni-Cox and multi-Cox regression analyses were performed to explore whether the prognostic signature was a potential independent prognostic indicator for patients with LUAD, and the results were visualized with two forest maps. Several receiver operating characteristic (ROC) curves were generated, and the area under the ROC curve (AUC) was calculated by the survival, survminer, and timeROC R packages to validate the predictive value of the prognostic signature.

We combined the risk score with the clinical variables of age, gender, N stage, T stage, M stage, and tumor stage to set up a nomogram for the 1-, 3-, and 5-year OS of LUAD patients by the rms R package. Correction curves based on the Hosmer–Lemeshow test were applied to illustrate the uniformity between the actual outcome and the signature prediction outcome.

We used gene set enrichment analyses (GSEA) software 4.1.2 (http://www.gsea-msigdb.org/gsea/index.jsp) to carry out GSEA and to identify significantly enriched pathways between the low- and high-risk groups. Values of p <0.05 and FDR <0.25 were considered the thresholds for statistical significance. The results were visualized by the gridExtra, grid, and ggplot2 R packages.

To figure out the relationship between this signature and TIME, firstly, we calculated the infiltration values for TCGA-LUAD dataset samples based on 7 algorithms: XCELL (Aran et al., 2017), TIMER (Li T et al., 2017; Li et al., 2020), QUANTISEQ (Finotello et al., 2019), MCPCOUNTER (Dienstmann et al., 2019), EPIC (Racle et al., 2017), CIBERSORT-ABS (Tamminga et al., 2020), and CIBERSORT (Chen et al., 2018). Using Spearman correlation analysis, the relationship of immune cell subpopulations and risk score value was evaluated. Wilcoxon signed-rank test, limma, scales, ggplot2, ggtext, tidyverse, and ggpubr R packages were applied, and the results are displayed in a bubble chart. Then, we explored the abundance of immune cells and stromal cells between different groups. The StromalScore, ImmuneScore, and ESTIMATEScore (StromalScore + ImmuneScore) of each patient were calculated. Their differences were compared using the Wilcoxon signed-rank test, and p <0.05 was considered to be significant. Subsequently, single-sample GSEA (ssGSEA) was conducted for scoring LUAD-infiltrating immune cells to quantify their relative content via the “GSVA” package. The scores of immune cells and pathways in different groups are shown on multi-boxplots, respectively. Finally, we also made comparisons about the immune checkpoint activation between low- and high-risk groups by the ggpubr R package.

The detailed flow diagram of our study is exhibited in Supplementary Figure S1. The transcriptome data of LUAD downloaded from TCGA included 59 normal samples and 539 tumor samples. We distinguished the mRNAs and lncRNAs by GTF files. According to the expression of 67 necroptosis genes and differentially expressed lncRNAs between normal and tumor samples, we finally obtained 1,016 NRlncRNAs (Supplementary Table S1), including 97 downregulated lncRNAs and 919 upregulated ones (Supplementary Figure S2).

Using uni-Cox regression analysis in the TCGA train set, we obtained 40 NRlncRNAs which were significantly correlated with OS and made a heat map (Figures 1A,B). To avoid overfitting and improve the accuracy and explainability of the prognostic signature, we performed the LASSO-penalized Cox analysis on these lncRNAs and extracted 19 lncRNAs related to necroptosis in LUAD when the first-rank value of Log(λ) was the minimum likelihood of deviance (Figures 1C,D). Subsequently, we constructed the predictive signature composed of 7 NRlncRNAs (AC026355.2, AC099850.3, AF131215.5, UST-AS2, ARHGAP26-AS1, FAM83A-AS1, and AC010999.2) via multi-Cox regression analysis. Of those lncRNAs, 6 lncRNAs were regulated positively by necroptosis genes in the Sankey diagram (Figure 1E). Meanwhile, some of those lncRNAs (AC099850.3, AF131215.5, and FAM83A-AS1) were demonstrated to be highly associated with NSCLC previously. Subsequently, the risk score of every LUAD patient was calculated based on correlation coefficients calculated by multivariate Cox regression analysis, and the patients were divided into low- and high-risk groups according to the median value of the risk score. The risk score was calculated as follows: risk score = (−0.3641 × AC026355.2 expression) + (0.1747 × AC099850.3 expression) + (−0.3943 × AF131215.5 expression) + (−0.6257 × UST-AS2 expression) + (−2.8454 × ARHGAP26-AS1 expression) + (0.3281 × FAM83A-AS1 expression) + (−2.1752 × AC010999.2 expression) (Supplementary Table S2).

To value the prognostic ability of the risk signature, we compared the distribution of risk score, the pattern of survival time, the survival status, and the relevant expression of 7 NRlncRNAs between the low- and high-risk groups in the train, test, and entire sets (Figures 2A–L). These all indicated that the low-risk group had better prognoses. Meanwhile, the LUAD patients were separated into groups according to age, gender, stage, T stage, N stage, and M stage to study the relationship between the risk signature and the prognosis of LUAD patients among universal clinicopathological variables. For different classifications, except T3-4 and M1 stage (Figures 3H, L), the OS of the patients in the low-risk group was significantly longer than that of the patients in the high-risk group (Figures 3A–G, Figures 3I–K). The possible explanation of the T3–T4 and M1 stage might be the limited number of patients due to poor prognoses in advanced NSCLC. These results suggest that the predictive signature can also predict the prognosis of LUAD patients in a different group of age, gender, stage, N stage, T1-2 stage, and M0 stage.

To determine whether the predictive signature is an independent prognostic factor for LUAD patients, Cox regression analysis was performed in the entire set. The Uni-Cox regression analysis showed that stage, T stage, N stage, and risk score were significantly associated with the OS of LUAD patients (Figure 4A). The multi-Cox regression analysis showed that only risk score (hazard ratio = 1.331, confidence interval = 1.175–1.507, p < 0.001) was an independent predictor of OS in LUAD patients (Figure 4B). Then, we used AUC to validate the sensitivity and the specificity of the signature in the entire set. The AUC of the risk score was 0.723, which was better than that of clinicopathological variables in predicting the prognosis of LUAD patients (Figure 4C). The AUCs of 1-, 3-, and 5-year survival were 0.723, 0.679, and 0.715, respectively, which indicated a good predictive value (Figure 4D). These results further implied that the signature was a promising biomarker for indicating the prognosis risk of LUAD.

The nomogram including clinicopathological variables and the risk score were constructed to predict the 1-, 3-, and 5-year prognosis of LUAD patients (Figure 5A). The calibration curves indicated a good consistency between the actual OS rates and the predicted survival rates at 1, 3, and 5 years (Figure 5B).

Based on the different prognoses of patients in the high- and low-risk groups, we conducted GSEA to explore the underlying differences in biological functions between risk groups. We found that the T/B cell receptor signaling pathway, Fc epsilon RI signaling pathway, cytokine receptor interaction, and JAK-STAT signaling pathway were significantly enriched in the low-risk group (Figure 6A), indicating that low-risk patients are closely related to tumor- and immune-related pathways. The GSEA results also revealed that the Notch signaling pathway, Wnt signaling pathway, and p53 signaling pathway, pathways in cancer and small cell lung cancer, were significantly enriched in the high-risk group. The Notch pathway plays a vital role in lung tumorigenesis and progression. Researchers found that cigarette smoke could promote LUAD progression via activating the Notch-1 pathway (Chiappara et al., 2022). Additionally, Notch-1 signaling synergized with Hif-1α could upregulate the expression of survivin in LUAD cell line A549 (Chen et al., 2011). The overexpression of Wnt pathway-activating genes and the down-expression of negative regulators of the pathway are closely correlated with NSCLC tumorigenesis, prognosis, and resistance to therapy (Stewart, 2014; Zeybek et al., 2022). The Wnt responder cells showed an increased tumor propagation ability, suggesting that they have features of normal tissue stem cells (Tammela et al., 2017). These mechanisms may explain why the high-risk group has a worse prognosis. Then, we studied the correlation between risk scores and tumor-infiltrating immune cells (Figure 6B). More immune cells are closely related to the low-risk group on different platforms. Consistently, we also found that StromalScore, ImmuneScore, and ESTIMATEScore in low-risk patients were significantly higher than those of high-risk patients (Figures 6C–E). To further explore the correlation between risk scores and immune cells and functions, we quantified the enrichment scores of ssGSEA for different immune cell subgroups, related functions, or pathways. The results exhibited that activated dendritic cells (aDCs), B cells, DCs, immature dendritic cells (iDCs), mast cells, neutrophils, T helper cells, T follicular helper (Tfh) cells, tumor-infiltrating lymphocyte (TIL), and T regulatory cells (Tregs) were significantly negatively correlated with the risk score (Figure 6F). Compared with the high-risk group, several immune pathways, e.g., checkpoint, cytolytic activity, human leukocyte antigen (HLA), T cell co-inhibition, T cell co-stimulation, and type II IFN response were higher in the low-risk group (Figure 6G). Furthermore, by comparing immune checkpoint activation between different risk groups, we found that almost all the immune checkpoints expressed more activity in the low-risk group, such as CTLA-4, HAVCR2 (TIM3), PDCD1 (PD-1), TIGIT, and CD70 (Figure 6H). These findings suggested that, in the low-risk group, the immunological function is more active and might be more sensitive to immunotherapy.