Identification and validation of a twelve immune infiltration-related lncRNA prognostic signature for bladder cancer

Identification of bladder cancer immune subgroups

A total of 405 bladder cancer patients from the TCGA-BLCA with 430 transcriptome profiles were obtained (Normal = 19, Tumor = 411). The transcriptome profiles of these samples were conducting ssGSEA analysis to evaluate the immune infiltration and immune related functions. Notably, there were 29 immunological marker gene sets, including immunological cell subtypes, pathways or functions associated with immune, have been used to perform the ssGSEA analysis. Then the samples were clustered into two groups according to the immune infiltration results, named accordingly, high immune cell infiltration cluster (Immunity H, n = 315) and low immune cell infiltration cluster (Immunity L, n = 96) (Figure 1A). The box chart also showed a significant higher Immune Score, Stromal Score, and ESTIMATE score in Immunity H group, while the Immunity L group showed a higher Tumor Purity (Figure 1B). Meanwhile, the expression values of CD274 (PD-L1) and HLA family was discovered significantly higher in Immunity H group than that in Immunity L group (Figure 1C and 1D, p < 0.001). Moreover, the CIBERSORT algorithm was used to verify the reliability of the unsupervised immune subgroups and found that there were higher infiltration of CD8⁺ T cells and M1 macrophages in the Immunity H group compared to the Immunity L group (Figure 1E). So, all these results showed that this bladder cancer unsupervised grouping performed well in distinguishing the difference of immune infiltration between samples and was suitable for further analysis.

Figure 1. Construction and verification of bladder cancer clustering by immune infiltration. (A) The immune cells were highly infiltrated in the cluster1 group, which was named as the Immunity_H (high immune cell infiltration group), and the low expression one in the cluster2 group was named as the Immunity_L (low immune cell infiltration group). The Tumor Purity, ESTIMATE Score, Immune Score and Stromal Score of each sample gene were also displayed with the grouping information by using ESTIMATE's algorithm. (B) The box-plot showed a statistical difference in Tumor Purity, ESTIMATE Score, Immune Score and Stromal Score between the two groups (p < .01). (C and D) In Immunity_H (red), the expression of HLA family genes and CD274 were all significantly higher than that in the Immunity_L (green) (p < 0.001). (E) The statistical chart showed the proportion difference of each immune cell between the Immunity_H (red) and the Immunity_L (green), after using the CIBERSORT method.

Identification of differentially expressed immune infiltration related lncRNAs

We firstly applied the criteria of |log₂FC| > 1 and FDR <0.05 to screen the differentially expressed lncRNAs between bladder cancer samples (n = 411) and paracancerous samples (n = 19). Thus, we discovered 1669 lncRNAs up regulated and 635 lncRNAs down regulated (Figure 2A). Later, the same criteria were conducted between the Immunity H/L groups, and 1601 differentially expressed lncRNAs were discovered (Figure 2B). Finally, we took an intersect of them and identified a total of 440 differentially expressed immune-infiltration-related lncRNAs (Figure 2C).

Figure 2. Analysis of the differentially expressed lncRNAs. (A) The volcano plot showed that 1699 and 635 genes were up-regulated and down-regulated between bladder cancer and paracancerous tissues. Each red dot showed an up-regulated gene and green showed downregulated genes (fold change >2, p < 0.05). (B) Consistent with Figure 3A, the volcano plot showed that 414 and 1187 genes were up-regulated and down-regulated between high and low immune cell infiltration group. (C) After taking an intersect, we obtained a total of 440 differentially expressed immune infiltration-related lncRNAs.

Construction and assessment of 12 immune-infiltration-related lncRNA prognostic signature for BC

We performed univariate Cox regression analysis of the 440 differentially expressed immune infiltration related lncRNAs. (Figure 3A). And 68 of them were found by applying the criterion of p < 0.05 (Original data sheet is in the Supplementary Table 1). Following this, we performed the least absolute shrinkage and selection operation (LASSO) regression to prevent overfitting and check 19 lncRNAs as appropriate variables (Figure 3B and 3C). Subsequently, we performed multivariate Cox regression to these 19 lncRNAs and developed a twelve-differentially expressed immune-infiltration-related lncRNAs signature (Figure 3D). The detailed information of these 12 lncRNAs was shown in Table 1 and the univariate regression results of these 12 lncRNAs were shown in Figure 3A. In addition, the risk score of each sample was calculated by the following formula: Risk score = 0.19 × AL136084.3 − 0.67 × AL590999.1 + 0.70 × AC090673.1 − 0.31 × AL078587.1 − 1.74 × AL096803.3 − 0.61 × AL357054.4 − 0.31 × AC073534.1 − 0.03 × PSORS1C3 − 0.15 × LINC02195 − 0.06 × AL731567.1 + 0.82 × AL022324.3 + 0.17 × AL591806.1 (Table 1). High-/low-risk groups were determined in samples according to the median risk-score. Kaplan-Meier curve showed that the overall survival (OS) of patients with high-risk is much lower than patients with low-risk, suggesting ab effective prognostic characteristics of risk score (p = 5.732e−11) (Figure 3E). The risk score distribution and corresponding scatterplot were together to show each bladder cancer sample’s risk and survival status. And we could see a higher risk score and mortality rate in high-risk groups (Figure 3F and 3G). Also, the heatmap suggested that there existed different expression between high-/low-risk groups (Figure 3H). In total, all these results confirm that this 12 immune-infiltration-related lncRNAs signature performs well in the prognosis prediction.

Figure 3. Identification and assessment of immune-related lncRNA prognostic signature for bladder cancer. (A) The HR and p-value of selected genes in the immune terms using the univariable Cox HR regression (Criteria: p-value <0.05). (B) The LASSO Cox analysis identified 19 lncRNAs most related to prognostics. (C) The 10-round cross-validation determined the optimal values of the penalty parameter. (D) The HR and p-value from the multivariable Cox HR regression prognostic signature. (E) Patients showed poor overall survival (OS) in the high-risk group (red) than those in the low-risk group (blue). (F) The risk curve of each sample reordered by risk score. (G) A sample survival overview using the scatter plot. The green dots represent survival and red represent death, respectively. (H) Heatmap showed the expression of the signature in the high-risk groups and low-risk groups. The pink and blue bars represented the low-risk group and the high-risk group. And the evolution from green to red represented the 0 to 6 level of gene expression.

12 immune-infiltration-related lncRNAs signature can be an independent prognostic factor in BC

Univariate and multivariate Cox regression were conducted to investigate if there was no association between 12 immune-related lncRNAs and clinicopathological factors or not. From the result, the HR (hazard ratio) and its 95% CI were 1.269 (1.206–1.336) in univariate Cox proportional regression (p < 0.001), and 1.244 (1.178–1.312) in multivariate Cox proportional regression (p < 0.001), separately, showing that this signature can perform as an independent prognostic element in bladder cancer patients (Figure 4A and 4B). Besides, time-dependent receiver operating characteristics (ROC) analysis was applied and plotted to compare the specificity and sensitivity of this signature with others clinical factors. The areas under the ROC curve (AUC) of the risk score in 1-. 3-, and 5-years were 0.741, 0.751 and 0.772, showing the 12 lncRNAs signature was more reliable than those commonly used clinicopathological factors including age, gender, grade, and stage (Figure 4C–4E). In total, both these two results suggested that the 12 lncRNAs signature could serve as an effective independent prognostic element for patients with bladder cancer.

Figure 4. Evaluate the independent prognostic value of the risk score by using the Cox regression analysis. The (A) univariate cox regression and (B) multivariate cox regression analysis of age, gender, grade, and risk score. Calculate the AUC for age, gender, grade, and risk score of the total survival according to the multivariate time-dependent ROC curve for 1-year (C), 3-years (D) and 5-years (E).

Relevance between 12 immune-infiltration-related lncRNAs signature and the immune cell infiltration

To investigate the correlation between this 12 immune-infiltration-related lncRNAs signature and the corresponding immune infiltration, we took another algorithm (TIMER) to assess the immune cell infiltration of each sample. Here in Figure 5A–5F, the infiltration of CD8+ T cells, dendritic cells and macrophages were significantly positive correlated with risk scores. Besides, we also used CIBESORT to analyzed the immune infiltration. And Figure 5G showed the more detailed relationship between different types of immune cells and risk scores, suggesting that the correlation co-efficient of macrophages M2 was also positive with risk scores, respectively, the same as the results showed in TIMER. In summary, these findings suggested that this 12-immune-infiltration related lncRNAs signature was associated with the immune infiltration.

Figure 5. Correlation between the 12 lncRNA prognostic signature for bladder cancer and the infiltration of immune cell subtypes. The correlation values of all the immune cells with risk score. (A) B cells. (B) CD4+ T cell. (C) CD8+ T cell. (D) Dendritic. (E) Neutrophil. (F) Macrophage. (G) The correlation of immune cell infiltration and risk scores according to the infiltration results estimated by CIBESORT.

Nomogram and drug response

It seems the risk score is the most weighted factor in the nomogram (Figure 6A), and all of the calibration curves show a consistency between the survival predicted by this nomogram and the actual survival rates (Figure 6B–6D). Besides, we performed drug sensitivity prediction. As shown in the Figure 7A–7E, the ordinate represents IC50, so the smaller the IC50, the more sensitive it is to drugs. So, it is interesting so to see that though patients with high risk scores are associated with a poor prognosis, they showed a more sensitive response to the cisplatin and doxorubicin than patients with low risk scores (Figure 7A and 7B), while they are less sensitive to the methotrexate (Figure 7C). The response to both gemcitabine and vinblastine are no differences (Figure 7D and 7E). Similarly, in our immunotherapy response prediction, the high-risk group showed significantly better response to anti-CTLA4 treatment, while the low-risk patients were more sensitive to PD1 monoclonal antibody (Figure 7F).

Figure 6. Quantitively analyze the overall survival of patients with bladder cancer. (A) Nomogram considering risk score and several common-used clinicopathological factors. Calibration curves for 1-yaer (B), 3-years (C), 5-years (D).

Figure 7. Prediction of the sensitivity of high risk patients to chemotherapy drugs and immunotherapy. (A) Cisplatin. (B) Doxorubicin. (C) Methotrexate. (D) Gemcitabine. (E) Vinblastine. (F) Immune checkpoint receptor: PD1-R, PD1-noR, CTLA4-R, ctal4-noR.

Data sources

TCGA_BLCA dataset were retrieved and downloaded in the fragments per kilobase of per million format (FPKM) from the TCGA database (https://portal.gdc.cancer.gov/), the corresponding clinical information containing age, gender, survival status, overall survival, pT stage, pN stage, pM stage, and AJCC stage were also downloaded from TCGA database. A total of 403 bladder cancer patients were enrolled in this research (Table 2). Notably, we obtained annotation gene sets consisting of 29 immune related gene sets considering both immune cells and immune related pathways or functions. Then we performed ssGSEA analysis to emphasize the integrative immune cells, immune related pathways, and immune related functions of each bladder cancer samples by the R package “GSVA”. Following this, all the bladder cancer samples were unsupervised clustered and divided into two clusters defined as high/low immune infiltration group according to the ssGSEA results.

Characteristics of the immune grouping

Having obtained the immune grouping of each sample, we would like to verify the effective ness of this immune cluster. Thus, we firstly conducted ESTIMATE analysis to emphasize the corresponding immune score, stromal score, estimate score, and tumor purity of each sample according to their transcriptional expression. Then, all these immune related scores were compared between these two immune groups. Besides, the gene expression of immune checkpoint including human leukocyte antigen (HLA) and CD274 (PD-L1) were compared between high immune cell infiltration group (Immunity H) and low immune cell infiltration group (Immunity L) to further verifying the effectiveness of the immune infiltration groups. Moreover, CIBERSORT algorithm was performed to estimate the detailed immune cell infiltration status of each sample, and Wilcoxon test was carried out to investigate the differential immune cell infiltration between these two immune infiltration groups.

Identification of differentially expressed immune-infiltration-related lncRNAs in BC

According to the immune infiltration groups clustered by the ssGSEA results, as mentioned above, the lncRNA expression profile data extracted and annotated from the transcriptome file were ranked from the Immunity H group to the Immunity L group. Then the immune infiltration related lncRNA were identified by the differential expression analysis between Immunity H group and Immunity L group. Besides, differentially expressed lncRNAs were identified by the same methodology between tumor tissue and normal adjacent tissue. Notably, both the filter criteria were the |log₂FC| >1 and FDR < 0.05. Following this, we took an intersection of these immune infiltration related lncRNAs and differentially expressed lncRNAs to obtain the final differentially expressed immune infiltration related lncRNAs for further analysis.

Further identification of immune infiltration-related lncRNA prognostic signature in BC

Having obtained the differentially expressed immune infiltration related lncRNAs, we merged the expression value with the detailed survival information of each sample. Then, we conducted univariate cox regression in order to screen if those differentially expressed lncRNAs have prognostic value or not. Besides, to avoid over-fitting, we subsequently carried the Least Absolute Shrinkage and Selection Operator (LASSO) regression to obtain the appropriated variables for further signature construction. Finally, multivariate cox regression was performed to construct the immune infiltration related lncRNA prognostic signature. A corresponding risk score formula was also established as follow:

$$riskScore={\displaystyle {\sum }_{i=1}^{n}coef(i)\cdot exp(i)}$$

Then each patient received a risk score according to this formula. The median value of all patients were set as the threshold, and all patients were divided into the high or low risk score group that higher than the median value is high risk and the lower represents low risk. Kaplan-Meier survival curves were plot and log-rank test was used to check whether this risk stratification is associated with the overall survival. Besides, univariate and multivariate cox regression were carried out to investigate whether this risk score could serve as an independent prognostic factor.

Correlation between prognostic signature and detailed immune infiltration

Here we took two differential methods to estimate the immune cell infiltration of each sample, and they were TIMER and CIBERORT. Among them, the immune cells calculated by TIMER was downloaded from the TIMER database (https://cistrome.shinyapps.io/timer/). The CIBERSORT results were calculated by the CIBERSORT algorithm. Following these two methods, PEARSON correlation test was carried out to further investigating the correlation between risk score and these detailed immune cells infiltration.

Nomogram and drug response

To emphasize the overall survival of patients with bladder cancer more quantitatively, we re-checked the corresponding clinical data of each patient and assembled a nomogram considering risk score and several common-used clinicopathological factors including age, gender, stage, and grade. Calibration curves for 1-, 3-, 5-years were also plotted to examine the accuracy of the nomogram. Finally, we separately predicted the drug response to chemotherapy by R package “ProPhetic” and the drug response to immunotherapy by submap algorithm [34]. Then compared the drug response to both chemotherapy and immunotherapy between high and low risk patients.

Statistical analysis

All statistical analysis was applied by R program version 4.0.2. Mean ± standard deviation was applied to describe the distribution of the continuous variables following normal distribution while the median (range) was used for continuous variables following abnormal distribution. Counts and percentages were used to describe the distribution of categorical variables. Two-tailed p < 0.05 was considered with statistical significance.

Data availability

Source data of this study were derived from the public repositories, as indicated in the section of “Materials and Methods” of the manuscript. And all data that support the findings of this study are available from the corresponding author upon reasonable request.