3.1 Identification of necroptosis-related lncRNAs in patients with BLCA
The study workflow is illustrated in Figure 1. In this study, we used data from 433 patients with BLCA from the TCGA cohort (T = 414, N = 19). After finding all lncRNAs and establishing coexpression analysis with 67 NRGs to identify lncRNAs related to necroptosis (coefficients >0.4 and p < 0.001), we established an interactive network to demonstrate the relationship between nrlncRNAs and NRGs (Figure 2A). Next, 686 nrlncRNAs with significant differences were identified by differential expression analysis (FC > 1 and FDR < 0.05). We found that 108 nrlncRNAs were downregulated, while the others were upregulated, mapping the volcano plot (Figure 2B). We selected up-and downregulated nrlncRNAs of the top 50 most significant differences in the heat map (Figure 3B).
3.2 Construction of prognostic necroptosis-related lncRNAs signature
According to univariate Cox regression analysis, nrlncRNAs related to survival (p < 0.05) were screened and plotted on forest map (Figure 3A). The Sankey diagram illustrates the regulatory relationship between these lncRNAs and NRGs (such as TSC1, BRAF, and ATRX) (Figure 3B). To avoid over-fitting of the signature, when the first value of Log (λ) was the minimum possible deviance, we performed LASSO regression and multiple Cox regression analysis on these nrlncRNAs (Figure 4A and B). Finally, an 11-necroptosis-related lncRNA signature was established.
Risk score was calculated using the following formula: Risk score = (AC008750.1 × -2.18520147058094)+（AL133297.1 × 2.05280592077303）+（AC005479.1 × 1.35382239745128）+ (AP003559.1 × 1.28533807015561) + (UBE2Q1-AS1 × -1.16093696050806) + (MIR100HG × 0.449749505941033) + (LINC02709 × 1.25314467417777) + (AC005387.1 × -0.607382896530388) + (LINC00649 × -0.490052688800663) + (AC021321.1 × -1.16608863558295) + (ETV7-AS1 × -0.589464760152836) (Table 1). Patients were further divided into high- and low-risk groups, with the median risk score as the standard. When the model was constructed, we divided the samples into two groups (train and test) and obtained OS analysis, risk score curve, survival status, and 11-nrlncRNAs risk heat map (Figure 5A–C).
Univariate and multivariate Cox regression analyses were used to compare the risk scores of the signature with age, gender, grade, and stage (Figure 6A and B). Area under the ROC curve (AUC) of 1-, 3-, and 5-years survival were 0.736, 0.735, and 0.744, respectively, and AUCs of the risk score (0.736) compared with other factors including age (0.662), gender(0.479), grade (0.529), and stage (0.641), respectively (Figure 6C and D).
We divided patients with BLCA into stage I/II and stage III/IV groups, and then compared whether OS (survival probability) of patients in high-or low-risk groups had significant differences in the two stages according to the risk score of the model. We built a nomogram for predicting the 1-, 3-, and 5-year survival rates of patients with BLCA and plotted the 1-, 3-, and 5-year calibration curves (Figure 6 E–I).
To investigate the functional and enrichment pathway differences between high-and low-risk groups, we used GSEA software to explore the Kyoto Encyclopedia of Genes and Genomes pathway throughout the entire set. The high-and low-risk groups enriched the first five pathways, and the BLCA pathway was built using multiple GSEA. The BLCA pathway (p 0.01 < 0.05; FDR 0.044< 0.25; |NES| 1.78> 1.5) was significant in the risk group of this model (Figure 7A and B) (Appendix Figure 1A)
3.5 Immunity factors, TME, and M6A analyses in risk groups
The relationship between immune cells and the patient’s risk score is depicted in the immune cell bubble chart (Figure 7C). Cancer-associated fibroblast_XCELL, endothelial cell_MCP-counter, granulocyte-monocyte progenitor XCELL, macrophage M1_QUANTISEQ, macrophage M2_CIBERSORT, etc. were positively correlated with BLCA patients’risk scores (p < 0.05). However, B cell memory CIBERSORT, B cell_TIMER, T cell CD4+ central memory_XCELL, T cell CD8+ naive_XCELL, and regulatory T cells (Tregs)_CIBERSORT-ABS, among others, showed a negative correlation(p < 0.05) (Figure 7C). Specific correlation coefficients and p-values for each type of immune cell are shown in Appendix Figure 1B and Appendix Table 3. SSGSEA showed that there were significant differences in immune cells such as macrophages, mast cells, neutrophils, Th1-cells, and Treg and immune-related function scores such as APC-co-inhibition and APC-co-stimulation between high- and low-risk groups, respectively (p < 0.001) (Figure 7D and E). Stromal cells, immune cells, and composite scores of each BLCA sample are shown, and the differences in these scores were compared between the high- and low-risk groups (p < 0.05) (Figure 7F–H). We found higher immune checkpoint gene expression in the high-risk group. There were significant differences in the expression of many genes, such as CD274, CD44, CD276, and HAVCR2, between the high-and low-risk groups (Figure 7I). Significant differences were found in multiple m6A regulators gene expression in risk groups such as FTO, ALKBH5, METTL3, and YTHDC1 (Figure 7G).
3.6 Distinguishing BLCA clusters and cluster evaluation
We performed nine cluster comparisons of BLCA samples and found the weakest correlation between the two subtypes. We divided the BLCA samples into two distinct subtypes (C1 and C2) (Appendix Figure 2). There are usually differences in immunotherapy and survival analysis between the C1and C2 subtypes (Figure 8A). We observed a significant difference in survival probability between C1 and C2 (p = 0.032) (Figure 8B). In addition, we used the ggalluvial graph to show the relationship between high/low risk groups and C 1/C 2 clusters (Figure 8C). PCA analysis showed that the expression of lncRNAs in BLCA could be used to scatter the sample clusters (C1 and C2) (Figure 8D and E).
The differences in stromal cells, immune cells, and composite scores in each BLCA sample were compared between the C1 and C2 groups (p < 0.05) (Figure 8F–H). Heat map revealed differences among different immune cells in the tumor clusters. T cell CD4+_TIMER, T cell CD8+_TIMER, and Neutrophil_TIMER immune cells were highly concentrated in the C1 cluster (Figure 8I). Immune checkpoint gene expressions such as TNFRSB4, CD274, and CD70 were mostly higher in the C1 cluster than in the C2 cluster (Figure 8J).
3.7 Drug sensitivity
The IC50 of many chemotherapeutic or targeted drugs showed significant differences between the high-and low-risk groups of the prognostic signature. This also suggests that the use of this model to differentiate patients can improve drug sensitivity to some extent and avoid drug resistance, especially some commonly used chemotherapeutic or targeted drugs such as cisplatin, dasatinib, docetaxel, gefitinib, imatinib, nilotinib, methotrexate, and lenalidomide (Figure 8A). The same pattern was predicted for the C1 and C2 clusters (Figure 8B).