Identification of immune-related LncRNA for predicting prognosis and immunotherapeutic response in bladder cancer

Data source and processing

As shown in Figure 1, the overlap was taken from 2420 IRlncRNAs and 1648 differentially expressed lncRNAs (DElncRNAs), and 190 differentially expressed IRlncRNAs (DEIRlncRNAs) were retained. Then, univariate Cox regression was conducted to choose characteristics for prognostic prediction of patients, and 42 DEIRlncRNAs with p<0.05 were retained for further analysis. The clinical characters of BLCA patients were downloaded from the UCSC database, and subsequently, TCGA dataset was split randomly into a training set (n=270) and a validating set (n=135) at a 2:1 ratio. There were no significant differences in age, gender, histologic grade, pathological stage, or diagnosis subtype between the two groups (Table 1). Then we identified 8 IRlncRNAs that were strongly associated with OS of BLCA by LASSO (Figure 2A, 2B) and multivariate Cox regression analysis in the training set. The forest plot shows the hazard radio (HR) with 95% confidence interval (CI) of eight IRlncRNAs (Figure 2C) and detailed information on these lncRNAs is listed in Table 2. Higher expression of MIR181A2HG (HR: 0.77 [95%CI 0.67-0.88], p<0.001), AC114730.3 (HR: 0.82 [95%CI 0.69-0.97], p=0.017) and LINC00892 (HR: 0.76 [95%CI 0.67-0.87], p<0.001) tended to predict increased survival, while higher expression of PTPRD-AS1 (HR: 1.13 [95%CI 1.01-1.27], p=0.036), LINC01013 (HR: 1.17 [95%CI 1.02-1.34], p=0.022), MRPL23-AS1 (HR: 1.12 [95%CI 1.03-1.22], p=0.007),LINC01395 (HR: 1.20 [95%CI 1.06-1.36], p=0.003) and AC002454.1 (HR: 1.13 [95%CI 1.00-1.28], p=0.046) tended to predict decreased survival.

An 8-IRlncRNA classifier to predict OS in BLCA

To assess the ability of the IRlncRNA signature to predict the survival of BLCA, we calculated the risk score for each case according to the expression of eight IRlncRNAs: Risk score = (-0.27 * expression value of MIR181A2HG) + (-0.20 * expression value of AC114730.3) + (-0.27 * expression value of LINC00892) + (0.12 * expression value of PTPRD-AS1) + (0.16 * expression value of LINC01013) + (0.12 * expression value of MRPL23-AS1) + (0.18 * expression value of LINC01395) + (0.13 * expression value of AC002454.1). Cases were split into high-risk and low-risk groups according to the median risk score (Figure 3A) and the mortality was higher in the high-risk group than in the low-risk group in the training set (Figure 3B). Moreover, MIR181A2HG, AC114730.3 and LINC00892 were highly expressed in the low-risk group, while PTPRD-AS1, LINC01013, MRPL23-AS1, LINC01395 and AC002454.1 were highly expressed in the high-risk group (Figure 3C). The results in the validating set were consistent with the findings described above (Supplementary Figure 1). Survival analysis revealed that patients in the high-risk group had a shorter OS than those in the low-risk group (p<0.001) in the training set (Figure 3D) and a similar result was observed in the validating set (p=0.002) (Figure 3E). Time-dependent receiver operating characteristic (ROC) curves were plotted and the area under curve (AUC) values of the classifier to predict 1-, 3- and 5-year overall survival were 0.72, 0.76, and 0.76, respectively (Figure 3F) in the training set and 0.74, 0.68, and 0.75, respectively, in the validating set (Figure 3G).

Survival prediction by the 8-IRlncRNA classifier is independent of clinical features

As shown in Table 3, clinicopathologic characteristics including age (p=0.038), and pathological stage (p=0.030) showed significant differences between the high-risk group and the low-risk group in the training set, while pathological stage (p=0.009) and diagnosis subtype (p=0.003) displayed distinct differences in the validating set. Patients with advanced pathological stage tended to obtain a high-risk score in the training (p=0.006) and validating groups (p=0.007) (Figure 4).

In the training set, the 8-IRlncRNA classifier, age and pathological stage were highly associated with OS by univariate and multivariate Cox regression analysis. Except for age, similar results were observed in the validating set (Table 4). Our results showed that 8-IRlncRNA classifier was an independent prognostic factor for OS in BLCA. A nomogram to predict 3- and 5-year overall survival utilizing the 8-IRlncRNA signature, pathological stage and histologic grade was developed (Figure 5A). The concordance index (C-index) of nomogram was 0.71, which increased the predictive power of OS compared with the 8-IRlncRNA classifier (C-index = 0.70). The calibration curves for 3- and 5-year overall survival showed that the predicted probability of OS was approximately equivalent to the actual OS (Figure 5B, 5C).

Prognostic value of the 8-IRlncRNA classifier for assessing recurrence-free survival

We also explored the prognostic value of the 8-IRlncRNA classifier to predict RFS of BLCA. Survival analysis showed that the RFS of patients in the high-risk group was significantly shorter than that in the low-risk group (Figure 6A, 6B). To evaluate the ability of the 8-IRlncRNA classifier to predict RFS of BLCA, we plotted ROC curves, and the AUC values of the classifier for prediction of 1-, 3- and 5-year RFS were 0.74, 0.7, and 0.71, respectively, in the training set (Figure 6C) and 0.67, 0.62, 0.68 respectively in the validating set (Figure 6D). Multivariate Cox regression analysis was conducted to identify prognostic factors for RFS, and outcomes indicated that pathological stage and the 8-IRlncRNA classifier were independent risk factors for RFS in BLCA patients (Supplementary Table 1). Our results demonstrated that this 8-IRlncRNA classifier could be a reliable prognostic predictor of BLCA recurrence.

Pathway enrichment analysis of the 8-IRlncRNA signature

To study the potential molecular mechanism of the 8-IRlncRNA classifier in BLCA, we performed gene set enrichment analysis (GSEA) of the high- and low-risk groups. The results showed that immune-related pathways such as antigen processing and presentation (Figure 7A) and hematopoietic cell lineage (Figure 7B) were inactivated in the high-risk group, suggesting that these pathways were correlated with the disease progression of BLCA. In addition, the top 5% of genes whose expression was correlated with the 8-IRlncRNA signature (p<0.05) were selected for enrichment analysis and some important immune-related pathways, including negative regulation of cytokine production, regulation of cytokine biosynthetic process, antigen processing and presentation, and PID CD8 TCR pathway were enriched (Figure 7C).

Correlation of the 8-IRlncRNA signature with immune cell infiltration

It was suggested that this 8-IRlncRNA classifier was related to immune-related pathways. Therefore, we explored the difference in immune cell infiltration between the two groups. Based on the ESTIMATE algorithm, we first calculated the stromal score and immune score of each BLCA sample. Higher stromal scores (-536.2 vs -816.8, p=0.005) and lower immune scores (-99.59 vs 151.2, p=0.009) were observed in the high-risk group compared with the low-risk group (Figure 8A, 8B), indicating different infiltration levels of immune cells in different risk groups. We further analyzed the abundance of 22 tumor-infiltrating immune cells in the two groups. In the high-risk group, the proportions of CD8 T cells (0.1105 vs 0.1371, p=0.014) and regulatory T cells (0.0207 vs 0.0384, p<0.001) were decreased, while the proportions of M0 macrophages (0.0764 vs 0.040, p=0.009), M2 macrophages (0.2323 vs 0.1880, p=0.002) and neutrophils (0.0075 vs 0.0031, p=0.009) were increased compared with those in the low-risk group (Figure 8C). Among these cells, low CD8 T cell infiltration was associated with low OS (Figure 8D, p=0.011) while high macrophage M2 cell infiltration was associated with low OS (Figure 8E, p=0.046). These findings strongly suggest that this IRlncRNA signature is associated with prognosis by interfering with immune cell infiltration in BLCA.

Potential of the IRlncRNA signature as an indicator of response to immunotherapy

Tumor immunotherapy using ICIs has been a promising treatment for advanced urothelial carcinoma [19]. It was confirmed that solid tumors deficient in MMR genes were usually immunogenic and showed extensively infiltrating T cells, making them highly responsive to ICIs [20]. We evaluated the correlation between the IRlncRNA signature and four key MMR genes (MSH6, MLH1, PMS2, MSH2). The results demonstrated that the risk score was significantly positively correlated with the expression of MLH1 (r=0.18, p<0.001) and MSH6 (r=0.23, p<0.001) (Figure 9A). The expression levels of MSH6 and MLH1 were both upregulated in the high-risk group (p<0.001 for MSH6, and p<0.001 for MLH1), indicating that this group benefited less from immunotherapy (Figure 9B). TMB was also a predictive biomarker for immunotherapy, and a high TMB suggested a high response rate to immunotherapy [21]. We acquired single nucleotide variation (SNV) data of 412 BLCA samples from TCGA and then selected the data processed by VarScan software for subsequent analysis. The landscape of mutation data in BLCA is showed in Supplementary Figure 2. We calculated the TMB for BLCA patients and matched the data for patients in our cohort. Patients in the high-risk group had a lower TMB than those in the low-risk group (4.313 vs 5.235, p=0.039), suggesting that these patients may be insensitive to ICIs (Figure 9E). Then, the association between the IRlncRNA signature and the expression levels of immune checkpoint genes (PD-1 and PD-L1) was investigated. As shown in Figure 9C, the risk score was significantly negatively correlated with the expression of PD-1 (r=0.17, P<0.001) and PD-L1 (r=0.352, p<0.001). The expression levels of PD-1 and PD-L1 were both downregulated in the high-risk group (p=0.003 for PD-1, and p<0.001 for PD-L1) (Figure 9D). These observed associations between our IRlncRNA signature and immunotherapy-related biomarkers indicated that BLCA patients in the high-risk group may be insensitive to ICIs. Therefore, the predictive value of the IRlncRNA signature was tested in the clear cell renal cell carcinoma (ccRCC) immunotherapy dataset [22]. The results demonstrated that patients in the nonresponse group had a higher risk score as defined by the 8-IRlncRNA signature (Figure 9F, p=0.046).

Data acquisition

Gene expression quantification data (FPKM and counts format) and SNV data for BLCA were downloaded from TCGA (https://portal.gdc.cancer.gov/). Then 19 normal samples and 411 BLCA samples were obtained. The matrix of RNA expression was extracted separately by annotations using the Gencode (GENCODE v 26) GTF file and normalized. Genes whose expression was “zero” in 90% of BLCA patients were removed. Clinical data were downloaded from the UCSC Xena website (https://xena.ucsc.edu/). To analyze the correlation of lncRNA expression signatures with the prognosis of bladder cancer patients, we filtered out samples without survival information. Then, we selected a total of 405 patients, and these patients were randomly divided into training (n=270) and validating sets (n=135) randomly at a 2:1 ratio for further analysis. Significant lncRNA-pathway pairs across 33 cancer types with each lncRNA having an activity in immune pathways (lncRES) score> 0.995 and a false discovery rate (FDR) < 0.05 were downloaded from Immlnc (http://biobigdata.hrbmu.edu.cn/ImmLnc/index.jsp) [38]. The list of immune-related lncRNAs in BLCA was extracted separately. Stromal scores and immune scores of BLCA were calculated by applying the ESTIMATE algorithm and downloaded from the website (https://bioinformatics.mdanderson.org/estimate/index.html) [39]. TMB was defined as the total number of somatic mutations per million bases and analyzed by the R package ‘maftools’. Because the ICIs treatment dataset was not available in BLCA samples, a dataset of ccRCC with available ICIs treatment data and transcriptomic profiles was obtained from the study of Miao et al [22]. Fifteen patients treated with PD-1 blockade therapy remained after excluding patients without expression of the 8 IRlncRNAs.

Analysis of differentially expressed lncRNAs

We obtained DElncRNAs between normal and tumor tissues, where P value <0.05 and |log2-fold change (FC)|> 1 were used as the cutoffs by using the R package ‘edgeR’ [40]. Then, we filtered DEIRlncRNAs by matching the list of immune-related lncRNA in BLCA. The R package ‘heatmap’ was used to display the eight selected IRlncRNAs.

Data processing and risk score calculation

First, the DEIRlncRNAs were subjected to univariable Cox regression analysis to select IRlncRNAs that were associated with the OS of BLCA patients. Final 42 IRlncRNAs with a p value<0.05 were found. Second, we conducted LASSO regression analysis to identify more meaningful prognostic variables. Finally, we used these IRlncRNAs in multivariable Cox regression to obtain the coefficients. Eight IRlncRNAs significantly correlated with OS were identified to build the prediction model weighted by their coefficients. A risk-score formula for OS was constructed, and each patient was assigned a risk score by this risk-score formula that was a linear combination of the expression levels of significant IRlncRNAs weighted by their respective Cox regression coefficients. The patients were divided into a low-risk group and a high-risk group based on the median risk score.

Identifying survival-related immune cells

Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSOFT) is an analytical tool utilizing deconvolution algorithm that can infer 22 human immune cell types and quantify the relative ratio of each cell type (http://cibersort.stanford.edu). To enhance the robustness of the results, CIBERSORT produces an empirical P value for the deconvolution of each sample based on Monte Carlo sampling [41]. RNA-Seq (FPKM format) data of BLCA was analyzed by R software to obtain the abundance ratio matrix of 22 immune cells in each sample. Samples with p > 0.05 were filtered out to increase the accuracy of the estimated results.

Pathway enrichment analysis

To explore the potential functions of the eight-IRlncRNA signature, we conducted GSEA to assess whether a predefined set of genes showed statistically significant, concordant differences according to the risk group by GSEA software (downloaded from http://software.broadinstitute.org/gsea/index.jsp) [42]. The gene database of “c2.cp.kegg.v6.2.symbols.gmt” from the molecular signature database was analyzed. Pathway enrichment analysis was conducted using the online database “Metascape” (http://metascape.org/) [43]. The significance threshold of FDR for enriched biological processes and pathways was set at 0.05.

Statistical analysis

The chi-square test or Fisher's exact test was conducted to measure the difference between the training and validating sets and the relationship between clinical data and risk score. Spearman's correlation coefficients were computed to investigate the potential relationship between two groups. Both univariable and multivariable Cox regression analyses were performed using the R package ‘survival’. The Kaplan-Meier survival curve with log-rank test was drawn to demonstrate the relationship between IRlncRNAs and OS or RFS by the R package ‘survival’. The Wilcoxon rank-sum test is a nonparametric statistical test mainly utilized for comparing two groups. The ROC curve was generated to measure the accuracy of survival prediction by the R package ‘survivalROC’. All statistical tests were two-sided, and p < 0.05 was considered statistically significant. All analyses were performed in SPSS version 25.0 (SPSS Inc., Chicago, IL, USA) or R version 3.5.2 (http://www.r-project.org/).