Constructing an immune‐ and ferroptosis‐related lncRNA signature to predict the immune landscape of human bladder cancer

Abstract Background LncRNAs play a variety of roles in the tumor microenvironment and cancer immune responses. Determining the significance of bladder cancer (BLCA)‐related genes to predict the prognostic and therapeutic response of BLCA is important. Methods IrlncRNA/ frlncRNA pairs were determined using univariate analysis. The signature was constructed based on this pairs. Finally, analysis and internal validation were performed from several aspects. Results We identified 60 immune‐ and ferroptosis‐related lncRNA pairs, among which 12 were included in the Cox proportional hazards model. Patients in low‐risk group survived for significantly longer. Survival and riskScore analyses showed that the low‐risk group had a significantly better clinical outcome. ROC curve analysis showed that AUC of OS values were more than 0.75 in the training set and the whole cohort. As assessed using Cox analysis, the riskScore was an independent prognostic predictor in the training, testing set and the whole cohort. The areas under the multi‐index ROC in the training set, the testing set, and the whole cohort were 0.777, 0.692, and 0.748, respectively. High‐risk group was positively associated with most of tumor‐infiltrating immune cells. High‐risk Scores correlated positively with high expression of CD274, but not with PD‐1. Low riskScores correlated positively with high expression levels of the genes ERBB2 and nectin‐4. High‐risk Score was associated with a lower IC50 value for Docetaxel, cisplatin, and Pazopanib, while there was an opposite result for metformin. Conclusions The signature constructed by pairing irlncRNAs and frlncRNAs showed a notable clinical predictive value.


| INTRODUC TI ON
Bladder cancer (BLCA) is the 10th most commonly diagnosed cancer worldwide, with approximately 573,000 new cases and 213,000 deaths in 2020 and is also the fourth most prevalent cancer in men. 1,2 Nonmuscle-invasive bladder cancer (NMIBC) and a part of muscle-invasive bladder cancer (MIBC) comprise non-metastatic bladder cancer. About a quarter of patients with BLCA suffer from MIBC or metastatic foci. 3,4 Moreover, the recurrence rate of BLCA is high, and after surgery, approximately 50% of patients suffer relapse and develop metastases. 5,6 Patients with BLCA have benefitted markedly from adjuvant chemotherapy and new immune checkpoint inhibitors (ICIs). [7][8][9][10][11][12] However, there is still a considerable proportion of patients who do not respond to immunotherapy at all stages of BLCA because of tumor immune evasion. 13 Lipid peroxidation-mediated and iron-dependent cell death is termed ferroptosis, which differs from autophagy, necrosis, and apoptosis. 14 Ferroptosis is involved in many diseases, including cancer. 15,16 In a recent study, researchers found that CD8+ T cells participate in the regulation of tumor ferroptosis during cancer immunotherapy. 17 Another study showed that immunotherapypromoted tumor ferroptosis and the effects of combination immunotherapy were enhanced using iron oxide-loaded nanovaccines (IONVs). 18 For bladder cancer, recently studies indicate that PGE2 metabolism affects variety of immune cells function and eventually lead to immune evasion. 19 In addition, the release of PGE2 could induct ferroptosis in cancer cells. 20 Altogether, these data indicate that ferroptosis is closely related to the antitumor immunity of bladder cancer.
Long non-coding RNAs (lncRNAs) are non-coding transcripts with a length >200 nucleotides. 21 LncRNAs comprise almost 80% of the human transcriptome and play a key role in post-transcriptional regulatory processes related to mRNA translation, stability, or splicing. 22 A recent study also showed that lncRNAs are critical to regulate genes encoding proteins that participate in cancer immunity. 23 LncRNAs are also important for immune-cell infiltration into tumors. 24 Previous studies have focused on signatures of immune-or ferroptosis-related lncRNAs, or on immune-and ferroptosis-related mRNAs. [25][26][27][28] Although these studies demonstrated excellent value for the prediction and prognosis for cancer diagnosis, evaluation, and treatment, the roles of immune-and ferroptosis-related ln-cRNA pairs are rarely studied. Therefore, the present study aimed to employ a novel algorithm to develop an immune-related lncRNA (irlncRNA) and ferroptosis-related lncRNA (frlncRNA) signature that does not rely on specific lncRNA expression levels and expected that it could more accurately predict the prognosis of patients and the response to immunotherapy in patients with bladder cancer. Its predictive value was estimated for patients with BLCA, and its chemotherapy efficacy, diagnostic effectiveness, immune checkpointrelated genes, and antibody-drug conjugate-related genes were determined.

| Transcriptome data retrieval, analysis of differential expression, and intersection analysis
We downloaded the human BLCA transcriptome profile (RNA sequencing data), harmonized as fragments per kilobase of transcript per million mapped reads (FPKM), from The Cancer Genome Atlas (TCGA: https://tcga-data.nci.nih.gov/tcga/) for follow-up analysis. To distinguish the lncRNAs from mRNAs, we also downloaded gene transfer format (GTF) files from Ensembl (http://asia.ensem bl.org). Then, we downloaded a list containing recognized immune-related genes (irgenes) from the ImmPort database (http://www.immpo rt.org), which was used in co-expression analysis to identify irlncRNAs. Similarly, coexpression analysis using a list of 288 recognized ferroptosis-related genes (fr-genes), downloaded from the FerrDb dataset (http://www. zhoun an.org/ferrd b/), was used to identify frlncRNAs. A correlation coefficient of more than 0.4 and a p-value <0.01 were used as the criteria to identify irlncRNAs and frlncRNAs. To identify the differentially expressed DEirlncRNAs and DEfrlncRNAs, the R package limma was used, with thresholds of log fold change (FC) >2 and a false discovery rate (FDR) <0.05. Venny 2.1 (https://bioin fogp.cnb.csic.es/tools/ venny/ index.html) was used to complete the gene intersection analysis.

| DE-IFRLs pairing
The DE-IFRLs were cyclically matched, and a 0-or-1 matrix was constructed supposing that C was a DE-IFRLs pair composed of lncRNA A and lncRNA B; If lncRNA A has a higher expression than lncRNA B, C was defined as 1; otherwise, C was defined as 0. Then, we further screened the established 0-or-1 matrix. The relationship between pairs was not considered if the expression level of an lncRNA pair was 0 or 1. The reason was that pairs could not predict patient survival outcome properly unless they had a certain rank. It was considered an effective matching if the number of lncRNA pairs with an expression level of 0 or 1 accounted for more than 20% of the total pairs.

| Patients' clinical data
The clinical data for patients with BLCA was retrieved from the BLCA project of the TCGA. Effective data were obtained by excluding repeated data and data with a follow-up of less than 30 days.

| Developing a risk model to evaluate the riskscore
First, we performed single factor analysis. Second, the patients were randomly divided into the training and testing sets at a ratio of 3:2 (237, 158) and we used least absolute shrinkage and selection operator (LASSO) regression in the training set, with a p-value of 0.05 and 10-fold cross validation. 1000 cycles of LASSO regression were run, and 1000 simulations were set for each cycle. Third, we recorded the frequency of each pair in the LASSO regression model. Fourth, Cox proportional hazard regression analysis was performed for pairs with a frequency >100, which were also used to construct the model. The following formula was then used to evaluate the riskScore for the constructed risk model for all the clinical cases: Clef (i) and E(i) represent the regression coefficient of the multivariate Cox analysis for the DEfirlncRNA pairs and the expression value of each DEfirlncRNA pair, respectively. Finally, patients in the training set were classified into low-and high-risk groups according to the median riskScore. The testing set and the whole cohort were also grouped using the same the median riskScore.

| Risk model validation
Univariate and multivariate Cox analysis were utilized to validate the lncRNA pairs in the model. Then, Kaplan-Meier analysis, time-dependent receiver-operating characteristic (ROC) analysis, riskScore analysis, and survival outcome analysis was used to demonstrate the survival difference of patients in the low-or high-risk groups. Next, R tools were used to visualize the results above by using the R packages of survival, survminer, and survivalROC.
Finally, to verify that the model could be used as an independent clinical prognostic predictor, multi-index ROC analysis, and univariate and multivariate Cox regression analyses were performed between the riskScore and the clinicopathological characteristics. The results are presented using a forest plot and ROC curve. These analyses used the R package survival.

| Tumor-infiltrating immune cell, immune checkpoint-related genes, and antibody-drug conjugate-related genes
To analyze the relationship between the riskScore and immune-cell characteristics, the immune infiltration status of the BLCA samples from the TCGA project were evaluated using CIBERSORT, 29,30 MCPcounter, 31 QUANTISEQ, 32,33 EPIC, 34 TIMER, 35,36 XCELL, 37,38 and CIBERSORT-ABS. 39 The relationship between the riskScore values and the immune infiltrated cells was determined using Spearman correlation analysis and a lollipop diagram was used to display the result. Then, the Wilcoxon signed-rank test was used to analyze the differences in immune checkpoint-and antibody-drug conjugaterelated genes among the different groups; the results of which are shown in a box chart. A p-value <0.05 was set as the significance threshold. The R limma, scales, ggtext, ggplot2, and ggpubr packages were used to perform these procedures.

| Determining the model's significance in clinical treatment
The BLCA dataset was mined for the half-maximal inhibitory concentration (IC50) values of commonly used chemotherapeutic drugs, with the aim of assessing the clinical significance of the model for BLCA treatment. The commonly used drugs comprised cisplatin, pazopanib, and docetaxel. Recent studies have suggested the therapeutic efficacy of metformin in certain cancers; therefore, we also explored the relationship between metformin and the model. The Wilcoxon signed-rank test was used to compare the differences in IC50 among the groups. The analysis was carried out using pRRophetic and ggplot2 of R and the results are presented using box plots.

| Establishing IFRL pairs and the risk assessment model
An iteration loop and a 0-or-1 matrix were used to screen 55 IFRLs, among which 1220 valid IFRL pairs were identified. A single factor test extracted 60 IFRL pairs, followed by modified LASSO regression analysis ( Figure 3A Figure 3C,D).

| Clinical evaluation using the risk assessment model
According to Kaplan-Meier analysis, the patients in the low-risk  Figure 5D,5F), but only the riskScore was

F I G U R E 2 Identification of differentially expressed immune-related lncRNAs (DEirlncRNAs) and ferroptosisrelated lncRNAs (DEfrlncRNAs) from TCGA data and Ensembl-based annotation. (A) DEirlncRNAs shown on a volcano plot (B) DEfrlncRNAs shown on a volcano plot. (C) Intersection of the two sets of DElncRNAs
an independent prognostic predictor in the testing set ( Figure 5E).

| Analyses of tumor-infiltrating immune cells, immunosuppressive molecules, and ADC targets using the risk assessment model
Immune-related genes and lncRNAs are interrelated; therefore, we explored the association between the model and the tumor immune microenvironment. The "infiltration_estimation_for_tcga.
csv" data file was downloaded from Timer database. We then performed correlation analysis utilizing Spearman analysis. The results showed that the high-risk group was positively associated with tumor-infiltrating immune cells, including monocytes, fibroblasts, and macrophages, but negatively related to CD4+ T cells and CD8+ T cells ( Figure 6A). ICIs comprise very important treatments for patients with BLCA in clinical practice; therefore, we assessed whether the model was associated with ICI-related biomarkers. The results showed that high-risk Scores correlated positively with high expression levels of the gene encoding CD274 (also known as PD-L1) (p < 0.05, Figure 6B); however, there was no significant relation between the riskScores and the expression levels of PDCD1 (also known as PD-1) ( Figure 6C). Antibody-drug conjugates represent a class of emerging therapeutics; therefore, we also investigated whether the model associated with relevant biomarkers. The results showed that low riskScores correlated (p < 0.001, Figure 6D) and nectin-4 (p < 0.001, Figure 6E).

| Analysis of the correlation between the risk model and chemotherapeutics
We attempted to discover the association between the riskScore and the efficacy of commonly used chemotherapeutic and targeted drugs. The results showed that a high-risk Score was associated with a lower IC50 value for Docetaxel (p < 0.001), cisplatin (p < 0.001), and Pazopanib (p < 0.001), which suggested that the developed model could be used to predict chemotherapeutic drug sensitivity ( Figure 7A-C). In the case of metformin, the IC50 in the low-risk group was lower than that in the high-risk group (p < 0.001, Figure 7D).

| DISCUSS ION
Immunotherapy of a considerable proportion of patients with BLCA is limited by immune invasion. In fact, recent studies have shown that, in vitro, macrophages effectively engulf ferroptotic cancer cells, supporting the existence of "find me" and "eat me" signals. 40,41 The Finally, the biological functions of the lncRNAs making up the prognostic signature need to be explored in detail in bladder cancer.

| CON CLUS ION
In the present work, we constructed an irlncRNAs and frlncRNAs signature that was independent of the expression levels of lncRNAs.
The signature could be used for prognosis prediction in patients with BLCA and could facilitate decisions regarding whether a patient might respond to BLCA immunotherapy and ADCs targeting Nectin-4 and HER-2.

ACK N OWLED G M ENT
We thank The Cancer Genome Atlas (TCGA) database for providing the data used in this study.

CO N FLI C T O F I NTE R E S T S
The authors have declared that no competing interest exists.

AUTH O R CO NTR I B UTI O N S
Xing Li and Jianting Xu conceived the study; Libin Zhou and Tefei Lu performed the R language analysis; Lei Zhang and Yanjun Li contributed significantly to the analysis and manuscript preparation; Xing Li performed the data analyses and wrote the manuscript; Huimin Long and Min Yin helped perform the analysis and provided constructive discussions.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets presented in the study are included in the article, further inquiries can be directed to the corresponding authors.