A Necroptosis-Related lncRNA to Develop a Signature to Predict the Outcome, Immune Landscape, and Chemotherapeutic Responses in Bladder Urothelial Carcinoma

Objective Many studies have drawn their attention to the immunotherapy of bladder urothelial carcinoma in terms of immunologic mechanisms of human body. These include immunogenicity of the tumor cells and involvement of long non-coding RNA (lncRNA). We constructed a necroptosis-related long noncoding RNA (nrlncRNA) risk factor model to predict BLCA outcomes and calculate correlations with chemosensitivity and immune infiltration. Methods Transcriptomic data from BLCA specimens were accessed from The Cancer Genome Atlas, and nrlncRNAs were identified by performing co-expression analysis. Univariate analysis was performed to identify differentially expressed nrlncRNA pairs. We constructed least absolute contraction and selector operation regression models and drew receiver operating characteristic curves for 1-, 3-, and 5-year survival rates. Akaike information criterion (AIC) values for survival over 1 year were determined as cutoff values in high- and low-risk subgroups. We reassessed the differences between subgroups in terms of survival, clinicopathological characteristics, chemotherapy efficacy, tumor-infiltrating immune cells, and markers of immunosuppression. Results We identified a total of 260 necroptosis-related lncRNA pairs, of which we incorporated 13 into the prognostic model. Areas under the curve of 1-, 3-, and 5- year survival time were 0.763, 0.836, and 0.842, respectively. We confirmed the excellent predictive performance of the risk model. Based on AIC values, we confirmed that the high-risk group was susceptible to unfavorable outcomes. The risk scores correlated with survival were age, clinical stage, grade, and tumor node metastases. The risk model was an independent predictor and demonstrated higher predictive power. The risk model can also be utilized to determine immune cell infiltration status, expression levels of immune checkpoint genes, and the sensitivity to cisplatin, doxorubicin, and methotrexate. Conclusion We constructed a novel necroptosis-related signature that predicts BLCA outcomes and performs satisfactorily in the immune landscape and chemotherapeutic responses.


INTRODUCTION
Bladder urothelial carcinoma (BLCA) is the 9 th most prevalent malignancy worldwide (1,2). More than 199,000 people died from the disease in 2018, and more than 549,000 new cases were confirmed in 2018 (3). The number of BLCA events has been increasing worldwide over the past two decades, and the incidence of BLCA has been increasing yearly due to population aging and environmental pollution (4,5). Treatment options for BLCA include transurethral resection, radical cystectomy, radiotherapy, and chemotherapy; nevertheless, BLCA remains an aggressive neoplasm with a substantial incidence of recurrence, metastasis, rapid progression, and unfavorable outcomes (6)(7)(8). One study in the United States estimated 80,470 new diagnoses of BLCA cases and at least 1,767 deaths in 2019 (9). Once the tumor has progressed and metastasized, the combination of systemic chemotherapy and surgery tends to be ineffective (10,11). At least 30% respond to immunotherapies and immune checkpoint inhibitors (ICIs) (12). There is a need for research to construct reliable prognostic biomarkers through molecular profiling to identify prognostic markers and therapeutic targets for BLCA to improve outcomes.
Necroptosis is programmed cell death triggered by MLKL, RIP1, and RIP3 (13,14). Several lines of evidence suggest that necroptosis is involved in Parkinson's disease, infectious diseases, cancer, and other diseases (14,15). Researchers found that necroptosis is a critical factor influencing tumor metastasis and T cell death (16). Interestingly, necroptosis has been related to antitumor immunity (13). As a substitute mode of programmed cell death to control apoptosis resistance, necroptosis plays a role in antitumor immunity in cancer therapy (13). Studies showed that necroptosis is a crucial cellular response that regulates many tumors' onset, progression, and metastasis (17). Studies also found that necroptosis can serve as a biomarker in some diseases, particularly cancer (18,19). Necroptosis enhances cancer cell migration and invasion in pancreatic carcinoma via the production of CXCL5 (20). Nevertheless, the precise functions of necroptosis in BLCA and the molecular mechanisms remain undetermined.
Long non-coding RNAs (lncRNAs) are a family of RNAs with no capacity for protein-coding that are more than two hundred nucleotides in length (21). A body of evidence suggests that necroptosis-related lncRNAs (nrlncRNAs) influence tumor progression and metastasis by triggering immune system processes and immune responses (17,22). Studies showed that nrlncRNAs are associated with outcomes of various tumors (18,22,23). Bioinformatics analysis based on The Cancer Genome Atlas (TCGA) revealed that lncRNAs are associated with the progression of BLCA via immune-related pathways (24).
Differentially expressed nrlncRNAs may serve as prognostic indicators and drug targets in BLCA. This study utilized gene expression profiles of high-throughput sequencing data from TCGA to identify lncRNAs targeting necroptosis-related genes and to develop necroptosis-related prognostic signatures for patients with BLCA.

Datasets and Preprocessing
We queried TCGA data portal (https://tcga-data.nci.nih.gov/ tcga/) (level 3 data, FPKM value) to obtain RNA sequencing profiles of 414 BLCA and 19 normal bladder specimens. We set the workflow type to "HTSeq-FPKM" and the data type to "Gene Expression Quantification" in the dataset download. To conduct subsequent analysis, normalization of the expression profiles to transcripts kilobase million values was carried out, and all analyses were conducted using R (version 4.1.1). Gencode (version 26) GTF files were obtained through Ensembl (http:// asia.ensembl.org) for annotation and differentiation of lncRNAs and mRNAs (25). Sex, age, clinical stage, and survival rates were obtained from TCGA for clinical data after removing specimens with insufficient clinical data or a survival duration of 0 days. Finally, we included 19 normal bladder specimens and 408 BLCA specimens. Supplementary Table 1 displays the clinical characteristics.

Identification of Necroptosis-Related lncRNAs
We accessed the Gene Set Enrichment Analysis site (http://www. gsea-msigdb.org/gsea/index.jsp) and obtained the necroptosis gene set M24779, which included eight necroptosis genes. We combined prior reports on necroptosis and obtained 67 necroptosis-related genes, details of which are displayed in Supplementary Table 2. We calculated Pearson correlations between necroptosis-related genes and identified lncRNAs. Correlation coefficients > 0.5 and p < 0.001 were used to identify lncRNAs linked to necroptosis.

Identification of Differentially Expressed Necroptosis-Related lncRNAs
We acquired 67 nrlncRNAs and used the R language version 4.1.1 "limma" package to identify differentially expressed necroptosis-related lncRNAs (DEnrlncRNAs) between BLCA and normal bladder specimens. The screening conditions were |log fold-change| > 1.0 and p < 0.05 (26). The expression matrix of differentially expressed lncRNAs was visualized using the heatmap package.

Paired DEnrlncRNAs
DEnrlncRNAs were identified using several pairing cycles, assuming that C was equal to the sum of lncRNA A and lncRNA B. A 0 or 1 matrix was then created. If lncRNA A expression level was greater than the level of lncRNA B, C was defined as 1; otherwise, C was defined as 0. The 0-or-1 matrix  was then re-evaluated. The correlation between pairs and patient outcomes was not evaluated if the expression level of lncRNA pairs was 0 or 1 because no pair ranks could accurately anticipate patient survival outcomes. In cases where the number of lncRNA pairs with an expression level of 0 or 1 accounted for more than 20% of the total logarithm, it was considered a valid match; otherwise, re-pairing was required.

Development of a Necroptosis-Related Prognostic Risk Model
According to BLCA data obtained from TCGA, univariate Cox proportional hazards regression analysis was utilized to identify lncRNA pairs associated with outcomes from necroptosis-related lncRNAs (p < 0.05). We then performed least absolute contraction and selector operation (LASSO) regression with 10-fold crossvalidation and a p-value of 0.05 and ran 1,000 loops. For each cycle, 1000 random stimuli were set to prevent overfitting. We then selected the best-paired combination to obtain nrlncRNA pairs in constructing the Cox risk coefficient model. By creating Cox univariate and multivariate analysis models, the risk coefficient for every necroptosis-related lncRNA pair associated with the outcomes of patients with BLCA was determined, and the risk score for each tumor specimen was determined. The aggregated risk score for every BLCA sample was the sum of the expression levels of each necroptosis-related lncRNA pair in the sample multiplied by the risk factor. The formula is as follows: The survminer and survival tools in R software were used to visualize the findings of the Cox analysis.

Assessing the Predictive Power of Prognostic Risk Model
The area under the curve (AUC) was evaluated to ascertain the predictive capacity of the risk model for determining patient outcomes, and receiver operating characteristic (ROC) curves were produced using the survivalROC module in the R software, which included ROCs at 1, 3, and 5 years. To compute the Akaike information criterion (AIC) value at every point on the 1-year ROC curve to obtain threshold values that optimize the aggregate of specificity and sensitivity in separating low-risk from high-risk individuals, we conducted a Kaplan-Meier analysis to identify disparities in survival between individuals in the high-and lowrisk groups, which we demonstrated using survival curves to calculate this cutoff value.

Prognostic Risk Model Validation
We performed the chi-square test to examine the correlation between the model and clinical and pathological features to assess the clinical significance of the constructed model. The Wilcoxon signed-rank test was performed for these clinicopathological variables to examine the difference in riskScore between groups. The analysis findings were displayed using box plots. We conducted univariate and multivariate Cox regression analyses between clinicopathological parameters and riskScore to determine whether the model could serve as an independent outcome predictor. To present the results, we generated forest plots. Survival, Heatmap, and ggupbr were the R packages we used.

Analysis Between Immune Checkpoints and Risk Models
To investigate the relationship between the expression of immune checkpoint-related genes (TIM-3, PDL1, LAG3 PD1, GAL9, TIGIT, PD1LG2, and CTLA4) and the model, we compared high-and low-risk subgroups and visualized them using the ggstatsplot package and violin plots.

The Value of Risk Models in Clinical Management
The half-inhibitory concentration (IC 50 ) of frequently used chemotherapeutic medicines for BLCA was determined in the dataset to assess the model's utility in medical therapy. The Wilcoxon signed-rank test was performed to calculate differences in IC 50 between high-and low-risk groups. The data were presented using the R software packages ggplot2 and pRRophetic (33).

Statistical Analyses
This study used R software (version 4.1.1) for statistical analysis. Differences between the two subgroups were estimated using the Wilcoxon rank-sum test. All statistical tests were two-way when p < 0.05 indicated statistical significance. Figure 1 demonstrates this study's flowchart. The initial step was to obtain transcriptome data for BLCA from TCGA. Finally, we included 19 normal samples and 414 BLCA samples. In the second step, data were annotated in accordance with the Gene Transfer Format files from Ensembl. We combined prior reports on necroptosis and obtained 67 necroptosis-related genes, according to 67 necroptosis-related genes and DEnrlncRNAs between normal and tumor samples (|log fold-change| >1.0 and p < 0.05), we identified a total of 291 nrlncRNAs (Supplementary Table 3), of which 89 were classified as DEnrlncRNAs ( Figure 2A); 76 underwent upregulation, and 13 u n d e r w e n t d o w n r e g u l a t i o n ( F i g u r e 2 B a n d Supplementary Table 4).

Development of DEnrlncRNAs Pairs and a Risk Model
Using multiple rounds of matching of 89 DEnrlncRNAs, 3101 necroptosis-related lncRNA pairs were identified (Supplementary Table 5). Next, univariate Cox regression analysis was performed to extract 260 DEnrlncRNA pairs affecting survival (Supplementary Table 6). To create a risk model, LASSO regression analysis was conducted to identify 13 necroptosis-related lncRNA pairs ( Figures 2C, D). Then, univariate and multivariate Cox regression analyses were conducted on these 13 nrlncRNA pairs ( Figures 2E, F), and each nrlncRNA pair's risk coefficient was calculated ( Table 1).

Evaluating the Risk Model's Outcomes Predictive Capability
The 13 nrlncRNA pairs were used to construct the 1-, 3-, and 5year ROC curves of BLCA patients ( Figure 3A), and the 1-year AUC was computed as a maximum of 0.763 ( Figure 3B). The AUC values for three and five years were 0.836 and 0.842, respectively, demonstrating that this risk model can also be used to predict 3-and 5-year outcomes for BLCA. According to best fit, the threshold value for differentiating between high and low-risk groups of BLCA patients was 1.189 ( Figure 3C).

Clinical Assessment by the Risk Model
We classified the patients into high-and low-risk groups based on the threshold value. We assigned 204 patients to the low-risk and 195 to the high-risk subgroup ( Figure 4A). Figure 4B demonstrates the distribution of survival status in the subgroups. The high-risk subgroup had more deaths than the low-risk subgroup, and the differences in survival time between subgroups were compared. Figure 4C illustrates that, compared to high-risk patients, low-risk patients had better outcomes (p < 0.001). We then conducted several chi-square tests to determine the relationship between risk scores and clinical and pathological features and generated a heatmap using the Wilcoxon signed-rank test ( Figure 5A). Scatterplots showed that survival ( Figure 5B ( Figure 5H), and N stage ( Figure 5I) were significantly associated with risk scores. Gender ( Figure 5D) was not significantly associated with risk scores. Then, the univariate and multivariate Cox regression analyses were performed on risk scores and clinical correlation factors, and forest maps were drawn ( Figures 6A, B). In univariate Cox analysis, we found that T stage, N stage, and riskScore were significantly associated with outcomes ( Figure 6A). Multivariate Cox regression analysis revealed that only riskScore could be utilized as an independent predictor for BLCA ( Figure 6B). ROCs were used to compare the differences in 1-year survival prediction performance. This risk model had the greatest AUC value (0.763) ( Figure 6C), suggesting that it has an excellent capacity to predict outcomes.   Table 7). Data visualization was performed using R language software ( Figure 7). The high-risk subgroup positively correlated with tumor-infiltrating immune cells such as cancer-related fibroblasts, T cell CD8+, M2 macrophages, and macrophage and was negatively associated with T cell CD4+ and T cell follicular helper cells (Supplementary Figure 1).  ICIs are therapeutic agents for management of BLCA. We determined whether the risk model was associated with biomarkers related to ICIs and found that PDL1 (p < 0.001; Figure 8A), HAVCR2 (p < 0.001; Figure 8C), LAG3 (p < 0.001; Figure 8D), PDCD1LG2 (p < 0.001; Figure 8G), and TIGIT (p < 0.001; Figure 8H) expression levels were significantly elevated in highrisk patients. Expression levels of CTLA4 (p > 0.05; Figure 8B) and PD1 (p > 0.05; Figure 8F) were increased; however, the differences were not significant. GAL9 (p < 0.001; Figure 8E) expression was attenuated in high-risk patients. These genes could be used as therapeutic targets for BLCA.

Correlation Analysis Between Risk Models and Chemotherapy Drugs
In addition to ICIs, chemotherapy is the first-line treatment for individuals with advanced BLCA. We also explored the correlation between risk models and the efficacy of conventional cancer medicines in BLCA. We found that a higher risk score was related to a reduced IC 50 for chemotherapy drugs including cisplatin (p < 0.001; Figure 9A) and doxorubicin (p < 0.001; Figure 9B), while it was associated with a higher IC 50 for methotrexate (p < 0.01; Figure 9D), gemcitabine ( Figure 9C), and vinblastine ( Figure 9E). Values were not significantly different between high-and low-risk groups, suggesting that the risk model predicts chemosensitivity.

DISCUSSION
An imbalance between tumor cell death and growth causes tumor formation and progression (34). Excessive cell growth or prevention of natural cell death exacerbates cancer progression. Some investigators argue that immortal cell proliferation and cell death suppression are distinct characteristics of malignant tumors (35). Necroptosis is a recently discovered type of cell death with morphological characteristics similar to necrosis. By contrast, necrosis refers to passive death induced by external physicochemical stress (e.g., inflammation or infection) and is not modulated by signaling pathways, while necroptosis is governed by programmed cell death (36). Many studies linked necroptosis to cancer incidence, progression, and metastasis (37,38). Necroptosis is also a viable strategy for eliminating cancer cells (39). Many studies investigated the role of lncRNAs in tumor onset and progression. Abnormally expressed lncRNAs in malignancies can be used as markers for clinical diagnosis, predicting outcomes, and developing targeted therapies (34). In BLCA, lncRNAs have been linked to cancer immunology and the tumor microenvironment (35). The lncRNA urothelial carcinoma-associated 1, the most studied lncRNA in BLCA, participates in several processes in the development of BLCA and is responsible for BLCA resistance (36,37). Based on the literature, several immune-related lncRNA models have been developed (24,38). These signature models were created using quantification of necroptosis-related lncRNA expression. To the best of our knowledge, there are no studies exploring the relationship between necroptosis-related lncRNAs and BLCA outcomes and underlying molecular mechanisms. We identified nrlncRNA pairs and developed a robust and independent risk profile of nrlncRNAs to determine the relationship between the model and BLCA outcomes and the potential impact on the BLCA tumor microenvironment and its corresponding treatment responses.
Zhang et al. evaluated the expression levels of ten hypoxiarelated lncRNAs to establish a signature predicting survival in BLCA (39). However, there is currently no study of nrlncRNAs in BLCA. Because of the critical role of nrlncRNAs, we developed a risk model with 13 DEnrlncRNAs pairs. This novel model is clinically useful, can distinguish high-or low-risk cases, and determine outcomes. In addition, we determined the correlation between the risk factor score and various clinical markers for each BLCA sample and found that the risk factor score independently predicted outcomes. We constructed ROC curves for clinically relevant indicators and compared the 1year ROC curves in the same chart. We confirmed that the risk factor score was the best predictor of BLCA outcome at 1 year, suggesting the robustness of the risk model. Then, we evaluated each of the DEnrlncRNAs identified in our model and found that they participate in the malignant phenotype of different cancer types, while lncSNHG14 overexpression promoted breast cancer proliferation and accelerated cell cycle progression (40). Wu et al.  found that lncRNA GHRLOS can be a biomarker for colorectal cancer metastasis and outcomes (41). LncRNA AC092279.1 was reported in a thyroid cancer prognostic model (42). LncRNA AC008115.3 was reported in a head and neck squamous cell carcinoma prognostic model (43). LncRNA AC010186.3 was reported in an ovarian Cancer prediction signature (44). LncRNA AC011503.2 and AC007128.1 were reported once in a bladder cancer predictive signature (45,46). Regarding AL136531.1, AC024361.1, AL731567.1, AC018521.6, AC024361.1, and GAS6-DT, there are no reports yet. Immune checkpoints and immune cell infiltration in cancerous tissues are critical for enhancing or inhibiting cell growth, invasion, and migration; immunotherapy is a novel therapeutic approach for managing diseases like BLCA (29). To explore the relationship between nrlncRNA signatures and immune cell infiltration, we used XCELL, TIMER, QUANTISEQ, EPIC MCPCOUNTER, CIBERSORT, and CIBERSORT-ABS algorithms to compare the content of immune cells in different risk score groups and found that the high-risk group was positively correlated with tumorinfiltrating immune cells such as cancer-associated fibroblast, T cell CD8+, M2 macrophages, and macrophages. It was negatively associated with T cell CD4+ and T cell follicular helper. Joźẃicki reported that breast cancer patients with reduced CD4+ T cell infiltration had shorter overall survival, and CD4+ T cells are a critical prognostic indicator, consistent with our findings (47). These findings suggest that this risk model could be used to anticipate the response to immunotherapy in patients with BLCA.
We also performed a correlation analysis of immune checkpoint genes and risk models and found that the PDL1, HAVCR2, LAG3, PDCD1LG2, and TIGIT expression levels were substantially higher in high-risk patients, and expression levels of GAL9 lower decreased in high-risk patients. These genes could be used as therapeutic targets for BLCA.
Bladder cancer is a complex malignant tumor, and chemotherapy is one of its essential treatment options. The guidelines recommend neoadjuvant chemotherapy before radical cystectomy of patients with BLCA, and the survival benefit of patients is close to 5-10%; nevertheless, some patients still do not respond to chemotherapy (48). Therefore, identifying predictors can avoid missing the optimal timing of surgery and minimize the adverse effects of chemotherapy. Here, we correlated BLCA chemotherapeutic agents with risk models and found that high-risk patients were more responsive to cisplatin and doxorubicin than low-risk individuals. Conversely, low-risk subjects were more responsive to methotrexate than high-risk subjects. These findings suggest that the risk model may help predict the sensitivity to doxorubicin, methotrexate, and cisplatin in patients with BLCA.
Although we employed rigorous approaches and algorithms to construct the model, this study has a few limitations. This study lacks external data to confirm the robustness of our risk model because the existing public databases do not include valid external data for verification. In a subsequent study, we will gather more clinical data and expand the sample size.
Finally, our findings revealed that a new signature created by nrlncRNAs might predict outcomes of BLCA and describe the immune landscape and chemotherapeutic therapy.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
The manuscript written was completed by JH, ZL, RD and GW, Experiment performance was done by GQ and YX, Data collection was conducted by HN, GY and CT, All the authors reviewed the manuscript and discussed the results and edited the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The Hunan Natural Science Foundation (#2021JJ50069) provided funding and grant its approval for the present research.