An immune-related lncRNA model for predicting prognosis, immune landscape and chemotherapeutic response in bladder cancer

Long noncoding RNAs (lncRNAs) participate in cancer immunity. We characterized the clinical significance of an immune-related lncRNA model and evaluated its association with immune infiltrations and chemosensitivity in bladder cancer. Transcriptome data of bladder cancer specimens were employed from The Cancer Genome Atlas. Dysregulated immune-related lncRNAs were screened via Pearson correlation and differential expression analyses, followed by recognition of lncRNA pairs. Then, a LASSO regression model was constructed, and receiver operator characteristic curves of one-, three- and five-year survival were established. Akaike information criterion (AIC) value of one-year survival was determined as the cutoff of high- and low-risk subgroups. The differences in survival, clinical features, immune cell infiltrations and chemosensitivity were compared between subgroups. Totally, 90 immune-related lncRNA pairs were identified, 15 of which were screened for constructing the prognostic model. The area under the curves of one-, three- and five-year survival were 0.806, 0.825 and 0.828, confirming the favorable predictive performance of this model. According to the AIC value, we clustered patients into high- and low-risk subgroups. High-risk score indicated unfavorable outcomes. The risk model was related to survival status, age, stage and TNM. Compared with conventional clinicopathological characteristics, the risk model displayed higher predictive efficacy and served as an independent predictor. Also, it could well characterize immune cell infiltration landscape and predict immune checkpoint expression and sensitivity to cisplatin and methotrexate. Collectively, the model conducted by paring immune-related lncRNAs regardless of expressions exhibits a favorable efficacy in predicting prognosis, immune landscape and chemotherapeutic response in bladder cancer.


Results
Identifying dysregulated immune-related lncRNAs in bladder cancer. Figure 1 depicted the workflow of this study. Here, transcriptome profiles of bladder cancer and normal specimens were obtained from TCGA and the lncRNAs were extracted. Immune-related lncRNAs that were distinctly correlated to immune-related genes were selected according to correlation coefficient > 0.4 and p < 0.001. As a result, 724 immune-related lncRNAs were identified (Supplementary table 1). Their expressions were compared between bladder cancer and normal specimens. Our data showed that 14 immune-related lncRNAs displayed down-    (Fig. 2C, D). Through uniand multivariate cox regression analyses, these pairs displayed distinct associations with survival outcomes (Fig. 2E, F). According to the coefficients and expressions of lncRNA pairs (Table 1), RS was calculated for bladder cancer subjects.
Evaluating the predictive performance of this risk model for prognoses. The cutoff value that differentiated bladder cancer subjects into high-and low-risk subgroups was 1.074 according to the AIC of one-year survival (Fig. 3A). The AUC of one-year survival was 0.806. This indicated the favorable predictive efficacy of this risk model. Furthermore, we conducted the ROCs of three-and five-year survival. The AUCs of three-and five-year survival were 0.825 and 0.828, demonstrating that this risk model was also utilized for predicting three-and five-year clinical outcomes of bladder cancer (Fig. 3B). According to the cutoff value, we clustered patients into high-and low-risk subgroups (Fig. 3C). The distributions of survival status between subgroups were depicted in Fig. 3D. High-risk subgroup possessed more dead patients in comparison to low-risk subgroup. The differences in survival duration were compared between subgroups. In Fig. 3E, low-risk patients were predictive of favorable clinical outcomes compared to high-risk patients (p < 0.001). Figure 4A depicted the associations between clinical features and this risk model in bladder cancer. We found that the risk model was in relation to survival status (p < 0.001), age (p < 0.05), stage (p < 0.05), T (p < 0.05), N (p < 0.05) and M (p < 0.05) of bladder cancer patients. The differences in RS were compared among different subgroups of clinical features. As shown in our data, patients with dead status exhibited higher RS than those with alive status (p < 2.22e−16; Fig. 4B). In Fig. 4C, patients in stage III-IV had elevated RS compared to those with stage I-II. Furthermore, > 65 patients displayed increased RS than ≤ 65 subjects (p = 0.006; Fig. 4D). As depicted in Fig. 4E, subjects with T3-4 possessed higher RS than those with T1-2. Compared to patients with N0, increased RS was detected in those with N1-2 (Fig. 4F). Also, higher RS was found in patients with M1 or Mx than M0 (Fig. 4G). There was elevated RS in high grade than low grade specimens (p = 0.0045; Fig. 4H). Nevertheless, no significant difference in RS was found between female and male specimens (Fig. 4I). Hence, this risk model might be relation to bladder cancer progression and metastases. This risk model as an independent prognostic predictor. As depicted in univariate cox regression analyses, stage, T, N, and risk model were in relation to bladder cancer prognoses (Fig. 5A). This was indicative that above factors might impact patients' clinical outcomes. To evaluate the predictive independency, multivariate cox regression analyses were conducted. In Fig. 5B, this risk model might be independently predictive of patients' prognoses. ROCs were conducted for comparing their differences in predictive performance of one- Table 1. Regression coefficients of each factor in this prognostic immune-related lncRNA signature. HR hazard ratio, HR.95L 95% CI lower limit, HR.95H 95% CI upper limit.   between this risk model and chemosensitivity were evaluated in bladder cancer specimens. Our data showed that high-risk patients exhibited decreased IC50 values of cisplatin in comparison to low-risk subjects ( Fig. 8A; p = 0.043). This indicated that high-risk scores were predictive of higher sensitivity to cisplatin. In Fig. 8B, reduced IC50 values of methotrexate were found in low-risk specimens than high-risk specimens (p = 1.5e−08), demonstrating that low-risk scores were in relation to higher sensitivity to methotrexate. Furthermore, this study evaluated the differences in IC50 values of vinblastine (Fig. 8C), gemcitabine (Fig. 8D) and doxorubicin (Fig. 8E) between high-and low-risk specimens. Nevertheless, no significant differences were found.

Discussion
LncRNAs have been confirmed to be in relation to cancer immunity and tumor microenvironment in bladder cancer 19 . Several immune-related lncRNA models have been constructed according to published literature 16,20,21 . Nevertheless, these signature models are developed on the basis of expression quantifications of immune-related lncRNAs. Herein, this study recognized the immune-related lncRNA pairs and constructed a reliable and independent risk model through combining two lncRNAs, not adopting their expression levels 22 .
Here, we firstly screened immune-related lncRNAs by Pearson correlation analyses. Different from previous research, dysregulated immune-related lncRNAs were further identified by comparing their expressions between bladder cancer and normal specimens 16,23 . With cyclically single pairing methods with 0-or-1 matrices, www.nature.com/scientificreports/ we identified immune-related lncRNA pairs. Combining univariate regression analyses and LASSO analyses, we developed a risk model for bladder cancer. Not using the median RS as the cutoff value that differentiated bladder cancer subjects into high-and low-risk subgroups, the AIC value of one-year survival was determined as the optimal cut-off value 16,24,25 . Furthermore, this risk model possessed distinct associations with survival status, age, stage and TNM of bladder cancer. According to multivariate regression analyses, the risk model might independently predict bladder cancer patients' OS. In comparison to other clinical features, the risk model displayed the highest AUC of one-year OS, indicating that the model possessed the potential as a favorable predictor of bladder cancer. Moreover, this risk model was in relation to immune cell infiltrations and immune checkpoints. Intertumoral tumor-infiltrating immune cells may impact the responses to ICIs 26,27 . Here, by comprehensively utilizing XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS and CIBERSORT algorithms, we characterized the correlations between risk model and immune cell infiltrations. High-risk specimens possessed elevated infiltration levels of myeloid dendritic cells, B cells native, macrophages M0 and M2, neutrophils and T cells CD8 in bladder cancer. Also, GAL9 displayed elevated expression in low-than high-risk specimens while higher TIM-3 and PD1LG2 expressions were found in high-risk specimens than low-risk specimens. These data indicated that this risk model might be utilized for predicting immunotherapy response of bladder cancer. Bladder cancer represents a complex malignancy correlated to high morbidity and mortality risks if not treated optimally. Neoadjuvant chemotherapy has been recommended prior to radical cystectomy for bladder cancer. Although the survival benefit is nearly 5-10%, some subjects cannot respond to chemotherapy 28 . Thus, identifying predictors may distinctly reduce side effects and miss the optimal time for surgery. Here, our data suggested that high-risk patients exhibited higher sensitivity to cisplatin in comparison to low-risk individuals. Inversely, subjects with low-risk were more sensitive to methotrexate than those with high-risk. Above data indicated that this risk model might possess the potential to predict the sensitivity to cisplatin and methotrexate for bladder cancer.
Due to high abundance, lncRNAs have distinct biological functions 29,30 . Our methods identified dysregulated immune-related lncRNAs and established the optimal immune-related lncRNA pairs. Hence, pairs with high or www.nature.com/scientificreports/ low expressions only were tested not detecting expression levels of each lncRNA. Our risk signature possessed the superiority in clinical practice for distinguishing high-and low-risk patients. Due to the closely correlations to immune-related genes, the selected lncRNAs potentially participated in modulating the immune microenvironment of bladder cancer. However, there are several limitations in our study. Firstly, because the lncRNA expression profiles of bladder cancer patients with complete survival time are not available in public databases, this study did not have an external validation to evaluate the performance of the prognostic immune-related lncRNA model. Therefore, more independent bladder cancer cohorts should be utilized for validating the risk model in our future studies. Furthermore, the functions of these lncRNAs and their interactions with immunerelated genes will be confirmed based on in vitro and in vivo experiments. Collectively, this prognostic signature constructed by 15 immune-related lncRNA pairs served as an independent predictor and displayed the favorable performance in predicting prognoses of bladder cancer. Also, it had the potential to predict immune landscape and chemotherapeutic response for bladder cancer patients.  Identifying dysregulated immune-related lncRNAs. Differential expression analyses of the immunerelated lncRNAs between tumor and normal bladder specimens were screened utilizing limma package 32 . The lncRNAs with |log fold-change|> 1.5 and false discovery rate (FDR) < 0.05 were screened. Above lncRNAs were visualized by heatmap package.
Pairing dysregulated immune-related lncRNAs. Cyclically singly paring dysregulated immunerelated lncRNAs were screened. The 0-or-1 matrices were developed if α = lncRNA-1 + lncRNA-2. α = 1 when www.nature.com/scientificreports/ lncRNA-1 expression was > lncRNA-2, while α = 0 when lncRNA-1 expression was < lncRNA-2. If the expression of lncRNA pair was 0 or 1, we thought there were no associations between this pair and prognoses, since the pair that did not a certain rank cannot be correctly predictive of patients' prognoses. If the number of lncRNA pairs that expression was 0 or 1 occupied > 20% of entire pairs, this was an effective match.
Establishing a prognostic risk model. Prognoses analyses of dysregulated immune-related lncRNAs were carried out through univariate cox regression models. LncRNAs with p < 0.05 could impact survival outcomes of bladder cancer. These lncRNAs were put into Least Absolute Shrinkage and Selector Operation (LASSO) model via glmnet package 33 . Penalty parameter tuning was carried out through ten-fold cross-verification. This analysis was run lasting 1000 cycles. The frequency of every pairing in the 1000-times-repeated LASSO model was retained and pairing with frequency > 100 times was chosen for constructing this model. Afterwards, a multivariate Cox regression model was conducted for determining the risk score (RS) through the coefficients and expressions of candidate lncRNA pairs according to the following formula: RS = k i=1 βiSi , where β represented the coefficient of lncRNA pair i and S represented the expression of lncRNA pair i.
Evaluating the predictive efficacy of the prognostic risk model. Receiver operator characteristic (ROC) curves were depicted for assessing one-, three-and five-year overall survival (OS). By calculating the area under the curve (AUC), the predictive efficacy of the prognostic risk model was determined. The Akaike information criterion (AIC) value of each point for the one-year ROC curves was calculated for identifying the maximum inflection point, which was selected as the cut-off value for distinguishing patients into high and low risk subgroups. Survival status of each subgroup was visualized. Prognoses analyses of high and low risk subgroups were conducted through Kaplan-Meier curves. Differences in survival were determined with log-rank tests.
Clinical feature assessment of the risk model. Associations between RS and clinical features (survival status, age, gender, grade, stage, T, N and M) were evaluated via chi-square tests. Also, RS was compared among different subgroups according to these clinical features. Univariate cox regression analyses were conducted for screening which factors could impact patients' survival. Hazard ratio and p values were separately calculated. Utilizing multivariate cox regression analyses, indicators that were independently predictive of survival were determined. One-year ROC curves were plotted for comparing the predictive performance of risk model and other clinical features.