Identification of a ten-long noncoding RNA signature for predicting the survival and immune status of patients with bladder urothelial carcinoma based on the GEO database: a superior machine learning model

Determination of prognostic lncRNAs

Figure 1 visualizes the identification process. After subjecting the lncRNA expression data to univariable Cox regression analysis by BRB-Array Tools, we identified 49 lncRNAs that strongly correlated with DSS (P<0.01). These 49 lncRNAs and their details, such as regression coefficients, P values and hazard ratios (HRs), are recorded in Supplementary Table 2. RSF-VH was performed on the expression profile consisting of these 49 lncRNAs. With this method, 10 lncRNAs remained for the construction of a prognostic signature for survival prediction. The details of these 10 prognostic lncRNAs, including their P values, regression coefficients, and HRs, are recorded in Table 1. In addition, the out-of-bag (OOB) importance value for each lncRNA is displayed in Table 1 and Supplementary Table 1. We observed that LOC105375787 was the strongest risk factor, while TTTY19 acted as the strongest cancer suppressor. The importance measures the increase (or decrease) in prediction error for the forest ensemble when a variable is randomly permuted in the OOB samples [14, 15]. Supplementary Figure 1 illustrates how the expression value of these genes and patient mortality correlate in random survival trees.

Figure 1.

Workflow of the construction and validation of the signature. ROC: Receiver operating characteristic; GSEA: gene set enrichment analysis; GEO: Gene Expression Omnibus.

Inferentially, positive coefficients implied that higher expression levels of 8 genes including LOC105375787, CYTOR, URB1-AS1, C21orf91-OT1, CASC15, LOC101928433, FLJ45139 and HOTAIR predicted shorter survival. In contrast, negative coefficients implied that higher expression levels of genes including LINC00960 and TTTY19 might predict longer survival. After grouping the patients using the k-means clustering algorithm based on the expression level of each lncRNA, we calculated Kaplan-Meier estimates to display the prognosis related to each prognostic lncRNA.

The overexpression of CASC15, URB1-AS1, FLJ45139, LOC105375787, HOTAIR was related to significantly shortened survival times in patients (log-rank test <0.05). High TTTY19 expression was significantly associated with longer survival in BLCA patients (log-rank test<0.05). The detailed results are shown in Supplementary Figure 2. Additionally, we found that the expression levels of CASC15, FLJ45139, LOC101928433 and C21orf91-OT1 (4 out of 8) were significantly increased in BLCA tissues compared to adjacent tissues or normal mucosas. The expression level of LINC00960 (1 out of 2) was decreased in BLCA tissues compared with normal tissues. However, in contrast to our expectation, lower expression was observed in BLCA tissues. This may be caused by the small sample size of normal bladder mucosa, as only 9 normal bladder tissue samples were included in the analysis (Supplementary Figure 3).

Risk formula and prognosis

To further investigate the association of this 10-lncRNA signature with BLCA prognosis, a prognostic model was constructed as follows:

∑_i Coefficient (lncRNAi) × Expression (lncRNAi), which was

Risk score=

(4.4 × expression level of LOC105375787) +

(1.2 × expression level of CYTOR) + (1.2 × expression level of URB1-AS1) +

(3.2 × expression level of C21orf91-OT1) +

(0.62 × expression level of CASC15) +

(0.91 × expression level of LOC101928433) +

(1.3 × expression level of FLJ45139) +

(1.9 + expression level of HOTAIR) +

(-0.39 × expression level of LINC00960) +

(-4.1 × expression level of TTTY19).

Using this formula, we calculated the risk score for each patient (Supplementary Tables 6–9). Defining the median risk score as a cut-off, the patients in the three series were divided into a high-risk subgroup or a low-risk subgroup of the same sample size (N=112, N=66, and N=83 and N=112, N=65, and N=82). Two prognostic nomograms for BLCA patients were successfully constructed to provide a clinically applicable quantitative approach for individual DSS prediction based on the signature and clinical prognostic factors, such as age, sex, grade and T stage (Figure 2A, 2D). The calibration curves in Figure 2B, 2C show a narrow margin between the predicted 3- and 5-year DSS rates and the actual values. The calibration plots in Figure 2E, 2F show excellent agreement between the predicted 3- and 5-year DSS rates and the actual observations. The C-index of the molecular nomogram was 0.88 (95% confidence interval (CI)=0.86-0.91) in GSE32894, 0.77 (95% CI=0.72-0.82) in GSE32548, and 0.73 (95% CI=0.68-0.77) in GSE13507. The C-index of the prognostic model combining the molecular signature with clinical risk factors was as high as 0.94 (95% CI=0.0.93-0.96) in GSE32894, 0.78 (95% CI=0.71-0.85) in GSE32548, and 0.90 (95% CI=0.87-0.91) in GSE13507. The condition number (K value) of the 10-lncRNA signature was calculated in every dataset which confirmed that there was no collinearity among these 10 genes (Supplementary Table 12). These results implied a superior predictive ability of the model, whether based on the 10-lncRNA signature alone or the signature combined with clinical parameters (Figure 3).

Figure 2.

Nomograms to predict 3- or 5-year DSS in the GSE32894 dataset. Instructions: locate each characteristic on the corresponding variable axis and draw a vertical line upwards to the points axis to determine the specific point value. Repeat this process. Sum the total points value and locate it on the total points axis. Draw a vertical line down to the 3- or 5-year DSS to obtain the survival probability for a specific bladder cancer patient. (A), Nomogram for predicting 3- or 5-year DSS in GSE32894 based on the 10-lncRNA signature. (B), Calibration curve for the prediction of 3-year DSS by the nomogram in (A). (C), Calibration curve for the prediction of 5-year DSS by the nomogram in (A). (D), Nomogram for predicting 3- or 5-year DSS in GSE32894 based on the 10-lncRNA signature combined with clinical risk factors. (E), Calibration curve for the prediction of 3-year DSS by the nomogram in (D). (F), Calibration curve for the prediction of 5-year DSS by the nomogram in (D).

Figure 3.

C-indexes with 95% confidence intervals of the prognostic models.

As illustrated in Figure 4A, the Kaplan-Meier analysis of GSE32894 showed that patients with higher risk scores had shorter DSS times than those with lower risk scores (log-rank test P < 0.0001). To further confirm the predictive value of this prognostic model, Kaplan-Meier analyses were performed to validate our signature in two independent external series, GSE32548 and GSE13507. An identical conclusion was obtained in GSE32548, where patients with high risk scores had significantly shorter DSS times than patients in the low-risk group (log-rank test P=0.00051) (Figure 4B). Similarly, in GSE13507, shorter DSS times were observed in the high-risk subgroup (log-rank test P=0.018) (Figure 4C). Even using overall survival (OS) as the follow-up endpoint, the risk score still separated patients with different survival times in the GSE13507 validation series (log-rank test P=0.025) (Figure 4D).

Figure 4.

(A–D) Kaplan-Meier survival curves in subgroups stratified by the risk score. The tick marks on the Kaplan-Meier curves represent censored subjects. The differences between the two curves were assessed by the two-sided log-rank test.

Scatter plots and heatmaps in GSE32548 were used to show the relationships among the risk score, the expression levels of ten lncRNAs and death from bladder cancer. Figure 5A shows that patients with higher risk scores tended to have higher mortality. The expression patterns of BLCA-specific lncRNAs were visualized in the heatmaps shown in Figure 5B. Patients with higher risk scores were observed to be more likely to have higher expression levels of risk-related lncRNAs, such as LOC105375787, CYTOR, URB1-AS1, C21orf91-OT1, CASC15, LOC101928433, FLJ45139, and HOTAIR. In contrast, patients with lower risk scores were more likely to have higher expression levels of suppressor lncRNAs (LINC00960 and TTTY19).

Figure 5.

(A) Patients’ survival status and DSS time. (B) Heatmap of the lncRNA expression profiles. Rows represent lncRNA expression, and columns represent patients. The middle dividing lines represent the median lncRNA risk score cut-off point. The graduated colour, from blue to red, represents the risk score.

Independence of the 10-lncRNA signature in survival prediction from tumour stage and tumour histopathological grade

A good prognostic biomarker should be independent of prevailing prognostic clinicopathologic factors. The traditional risk stratification of survival mainly depends on histopathological evidence, such as tumour stage and histopathological grade. To confirm the independence and applicability of our 10-lncRNA signature, multivariate Cox regression along with stratification analysis were conducted in GSE32894 and GSE32548.

Figure 6A indicates that both the 10-lncRNA risk score (HR=1.55, 95% CI=1.223–2.0, P <0.001) and tumour stage (HR=12.04, 95% CI=2.665–54.4, P < 0.001) are predictors independent of age, sex and histopathological grade in GSE32894. In GSE32548, using age, sex, tumour stage, histopathological grade and carcinoma in situ (CIS) as the covariables, multivariable Cox regression analysis further confirmed that the 10-lncRNA risk score (HR=1.90, 95% CI=1.201–3.0, P <0.007) was a prognostic predictor independent of tumour stage and histopathological grade (Figure 6B).

Figure 6.

Forest plots showing the hazard ratios (HRs) with 95% confidence intervals (95% CIs) based on the multivariate Cox regression results. (A) Risk score and tumour stage are dependent of age, sex and histopathological grade in GSE32894. (B) Risk score is dependent of sex, age, tumour stage, histopathological grade and concomitant CIS in GSE32894.
(C) Risk score is dependent of sex, age and pN in patients with cystectomy in GSE32894. (D) Risk score is dependent on sex, age, histopathological grade, multifocality and tumour size in patients with NMIBC in GSE32548.

Figure 6.

Forest plots showing the hazard ratios (HRs) with 95% confidence intervals (95% CIs) based on the multivariate Cox regression results. (C) Risk score is dependent of sex, age and pN in patients with cystectomy in GSE32894. (D) Risk score is dependent on sex, age, histopathological grade, multifocality and tumour size in patients with NMIBC in GSE32548.

Stratification analysis was used to investigate whether the signature could discriminate patients with different prognoses irrespective of the same tumour stage. After stratifying the patients from GSE32894 and GSE32548 into the muscle invasion (≥T2) subgroup or non-muscle invasion (Ta or T1) subgroup, Kaplan-Meier plots showed that for patients in the same tumour invasion subgroup, those with higher risk scores had significantly shorter DSS times than those with lower risk scores (Figure 7A, 7B).

Figure 7.

Kaplan-Meier estimates along with stratification of DSS for patient grouping by the signature. (A) Kaplan-Meier curves for GSE32894 (N = 224), stratified by whether muscle invasion exists (≤T2 or >T2). (B) Kaplan-Meier curves for GSE32548 (N = 131), stratified by whether muscle invasion exists (≤T2 or >T2). (C) Kaplan-Meier curves for GSE32894 (N = 224), stratified by tumour grade. (D) Kaplan-Meier curves for GSE32548 (N = 131), stratified by tumour grade. (E) Kaplan-Meier curves for patients with cystectomy in GSE32894 (N = 51), stratified by pN (pN- or pN+). (F) Kaplan-Meier curves for patients with NMIUC (G) Kaplan-Meier curves for patients with NMIUC (3 cm). (H) Kaplan-Meier curves for GSE32548 (N = 131), stratified by concomitant CIS; the tick marks on the Kaplan-Meier curves represent death from disease. The differences between the two curves were assessed by the log-rank test. MI: muscle invasion; NMI: non-muscle invasion.

Considering that tumour grade is another primary prognostic factor in BLCA, stratification analysis was conducted to investigate whether the 10-lncRNA signature could accurately predict the prognosis of bladder cancer patients with the same histopathological grade. The patients in both GSE32894 and GSE32548 were classified into three groups (G1, G2 and G3) based on different tumour grades. Kaplan-Meier plots showed that for patients with the same tumour histopathological grade, patients with higher risk scores had significantly lower DSS curves than those with lower risk scores (Figure 7C, 7D).

Independence of the 10-lncRNA signature in survival prediction from lymph node metastasis, multifocality, tumour size and concomitant CIS

Other clinical and pathological characteristics, such as lymph node metastasis, multifocality, tumour size and concomitant CIS, have been widely recognized to be the prevailing predictors for the prognosis of BLCA patients. Multivariate Cox regression along with stratification analysis were then conducted to determine the prognostic independence of the signature.

Fifty-one patients in the GSE32894 series had available post-cystectomy information, which made it possible to investigate whether the signature is independent of lymph node metastasis. We performed multivariable Cox regression analysis on this cohort, including age, sex, pN at cystectomy, and the 10-lncRNA risk score as the covariables. The analysis showed that the 10-lncRNA risk score (HR=1.4, 95% CI=1.10–1.80, P < 0.007) was the only independent prognostic factor and had a close correlation with DSS (Figure 6C). Kaplan-Meier plots further confirmed that patients with higher risk scores had shorter DSS times than those with lower risk scores, despite having the same pN at cystectomy (Figure 7E).

Clinical information on multifocality and tumour size was available for 93 non-muscle-invasive urothelial carcinoma patients (tumour stage Figure 6D). Since the sample size of non-muscle-invasive urothelial carcinoma was small, Kaplan-Meier plots generated indistinctive Kaplan-Meier curves (P =0.065 and 0.49) among the stratifications comparing multifocality and tumour size. However, patients with higher risk scores were still more likely to have a poor prognosis, which was indicated by the red or green curves in Figure 7F, 7G.

When compared with concomitant CIS, age, sex, tumour stage and tumour grade, the 10-lncRNA risk score (HR=1.90, 95% CI=1.201–3.0, P <0.007) was determined to be a survival predictor independent of concomitant CIS according to previous results (Figure 6B). Kaplan-Meier plots with stratification showed that the 10-lncRNA signature was able to classify patients who had different survival times and the same concomitant CIS status (Figure 7H) at some level, though the P value did not reach the threshold (log-rank test P=0.06>0.05).

The prognostic value of the 10-lncRNA signature

To determine the predictive power of this 10-lncRNA signature, ROC curves were employed, and AUC values were calculated. AUCs were compared between the signature and other prevailing prognostic predictors.

In all three independent series, the 10-lncRNA signature demonstrated a discriminatory ability for predicting DSS, with AUC values of 0.871 (95% CI=0.808–0.934) in GSE32894 (Figure 8A), 0.752 (95% CI=0.624–0.879) in GSE32548 (Figure 8B) and 0.707 (95% CI=0.597–0.816) in GSE13507 (Figure 8C). As shown by our data, this 10-lncRNA signature had robust sensitivity and specificity to potentially become a superior prognostic biomarker.

Figure 8.

ROC curves showing the prognostic performance of the 10-lncRNA signature compared with that of different prognostic predictors. (A) Comparison of tumour stage and tumour grade in GSE32894. (B) Comparison of tumour stage and tumour grade in GSE32548. (C) Comparison of tumour stage and tumour grade in GSE13507. (D) Comparison of multifocality and tumour size in patients with NMIUC from GSE32548. (E) Compared with concomitant CIS in GSE32548. (F) Comparison with pN in patients with MIUC from GSE32894. NMIUC: non-muscle-invasive urothelial carcinoma; MIUC: muscle-invasive urothelial carcinoma.

In GSE32894, as shown in Figure 8A, both the 10-lncRNA risk score (AUC=0.871) and tumour stage (AUC=0.920) showed a high predictive performance with no significant difference (P=0.134). Compared with histopathological grade, the 10-lncRNA risk score showed greater prognostic potential (0.871 versus 0.772, P=0.0038). In GSE32548, as shown in Figure 8B, since the P value of DeLong's test did not reach the threshold, we cannot tell whether the risk score has a greater capability to predict survival than tumour stage (0.752 versus 0.824, P=0.306) or histopathological grade (0.752 versus 0.675, P=0.239). However, with an absolute AUC value estimated to be 0.752, the prognostic value of this 10-lncRNA signature was still considered remarkable. In GSE13507, as presented in Figure 8C, despite a smaller AUC (AUC=0.707), which was still estimated to exceed 0.70, there was no significant difference between the signature and tumour stage or histopathological grade (P=0.05 and P=0.36, respectively).

ROC curves were also calculated for the 10-lncRNA risk score, concomitant CIS, lymphovascular invasion, multifocality and tumour size. Compared with tumour multifocality and tumour size, as presented in Figure 8D, the AUC of the 10-lncRNA risk score was roughly the same (0.637 versus 0.750; 0.637 versus 0.685, 95% CI=0.393–0.882, P=0.5667; P=0.9366). As shown in Figure 8E, the AUC of the 10-lncRNA risk score was estimated to be greater than 0.70 in GSE32548, which was larger than that of concomitant CIS despite P>0.05 (0.752 versus 0.576, 95% CI=0.624–0.879, P=0.5822). Compared with lymphovascular invasion, as shown in Figure 8F, the 10-lncRNA risk score performed as a relatively good predictor and had a larger AUC (0.725 versus 0.613, 95% CI=0.574–0.876 P=0.22), albeit with no significance observed.

Relationship among the 10-lncRNA signature, potential therapeutic target signature and immune-checkpoint blockade (ICB) immunotherapy-related signature

A growing body of evidence highlights the key role of immune regulation in BLCA, involving some immune checkpoints, such as CTLA4, PD-1, and PD-L1. In addition, several genomic alterations, such as HER-2 (ERBB2), ERBB3 and FGFR3, which were included in this study, have been identified to be amenable in principle to therapeutic targeting and were reported to be associated with advanced BLCA in previous studies.

Here, we compared our 10-lncRNA signature with these biomarkers to determine their prognostic value and potential interconnection. Pearson correlation coefficients demonstrated that the 10-lncRNA signature had a negative correlation with FGFR3 and a positive correlation with CD274 (PD-L1) and CTLA4 (P<0.01) (Figure 9A). The AUC values of the ROC curves were calculated for each biomarker. The ROC curves showed that the 10-lncRNA signature had the largest AUC (estimated to be 0.871) compared with these therapeutic targets and immune checkpoints (DeLong's test P < 0.05) (Figure 9B). Assuming these molecules are all promoters regulating cancer development, the results indicated that our lncRNA signature had better stability and reliability in predicting the DSS of patients with BLCA and implied a closer relation between the cancer and the signature.

Figure 9.

ROC curves and correlation analysis of different prognostic signatures. (A) Correlation analyses among known immune checkpoints, potential therapeutic targets and the risk score. The circle size represents the P value, and the colour represents the correlation coefficient. Blue indicates a positive correlation, and red indicates a negative correlation. (B) ROC curves show the sensitivity and specificity of our 10-lncRNA signature and other known biomarkers in predicting the DSS of patients from GSE32894.

Determination of disease-related pathways

After grouping the patients in GSE32894 into two groups by the median risk score, the whole-genome expression profiles were subjected to analysis with GSEA 3.0 software to screen for pathways related to the signature. Figure 10 visualizes the enriched canonical pathways that were identified as “enriched (FDR < 0.01; P<0.05)” and maps them in a whole pathway overview. The immune system, DNA repair and cell cycle were observed to be the top three enriched pathways. Activation of E2F1 target genes at G1/S, G1/S-specific transcription, activation of the pre-replicative complex, activation of ATR in response to replication stress, and scavenging by class A receptors were the top 5 downregulated canonical pathways. RNA polymerase I promoter clearance, glucuronidation, RNA polymerase I transcription initiation, RNA polymerase I promoter escape, and synthesis of glycosylphosphatidylinositol (GPI) were the top 5 upregulated canonical pathways. Supplementary Table 10 contains all these significantly enriched canonical pathways.

Figure 10.

GSEA identified biological processes that correlated with the 10-lncRNA signature. Enriched canonical pathways (FDR < 0.01) were mapped in the whole pathway overview. Each dot represents a function, and each line indicates a pathway.

Analysis of the immune status of the high-risk and low-risk populations

To verify whether the 10-lncRNA signature can reflect the status of the tumour immune microenvironment, we analysed the relationships between the signature and immune cell infiltration. Figure 11A shows that the proportions of infiltrating B cells, naive CD4+ T cells, naive CD8+ T cells, dendritic cells, exhausted T cells, gamma delta T cells, macrophages, monocytes, neutrophils, NK cells, natural regulatory T cells, and follicular helper T cells were significantly enriched in the high-risk group (P<0.05). Central memory T cells, effector memory T cells, mucosal-associated invariant T cells, NKT cells, T helper 17 cells, and T helper 2 cells were significantly enriched in the low-risk group (P<0.05). The infiltration score retrieved from the ImmuCellAI resource is a reflection of the overall tumour immune microenvironment. Figure 11B shows that tumours in the low-risk group have higher infiltration scores, which suggests the connection between the prognostic signature and the inflammatory environment of bladder tumours (P<0.05).

Figure 11.

(A) The difference in tumour-infiltrating immune cells among risk groups as defined by the 10-lncRNA prognostic signature. (B) The difference in infiltration score among risk groups as defined by the 10-lncRNA prognostic signature.

Material preparation

We downloaded three whole-genome expression microarray data series from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). All data series were included based on the following criteria: 1. public expression data generated for BLCA were obtained; 2. the same manufacturer platform was used (in this paper, the following Illumina human expression beadchip platforms were used: GPL6947 and GPL6102); and 3. raw nonnormalized data and matched clinical data with follow-up information were obtained. Ultimately, GS32894, GSE32548 and GSE13507 were included in the present study after an initial quality check. The clinical information of the samples included in this study from GSE32894, GSE32548 and GSE13507 was recorded in Supplementary Tables 3–5. The information relative to the selected datasets has been organized in Supplementary Table 11. The BLCA samples in GSE32894 were defined as a training series (N=224), and the samples in GSE32548 were defined as an independent test series (N=131). In addition, 165 samples from GSE13507 were defined as an external validation series (N=165). In this study, we set two independent external datasets (GSE32548 and GSE13507) as verifiable groups rather than separating them because multiple independent validation datasets with larger sample sizes could reduce disease heterogeneity and improve model stability.

Data pre-processing and lncRNA mining

All bioinformatic approaches used in the present article were conducted with R software version 3.5.2 [50]. The raw nonnormalized expression data of the three series that were directly downloaded from the GEO database underwent noise correction and quantile normalization using the “limma” package in R. This method also log2 transformed the expression values and filtered control probes, leaving only the regular probes [51].

In GEO or other databases, we failed to find a microarray specifically designed for lncRNAs with the requisite clinical information. Therefore, by annotating the probes, we separated a lncRNA profile from total RNA expression data. The gene sequence IDs obtained from the GPL6947 platform mainly contained Unigene ID and Refseq_ID. The “org.Hs.eg.db” package in R was used to obtain three maps of the correspondence between diverse gene sequence IDs, including Unigene ID - Entrez ID, Refseq_ID - Entrez_ID and Entrez ID - Ensembl ID. According to these maps, each Unigene ID and Refseq_ID was mapped to the corresponding Entrez ID. After that, the Entrez ID was further matched to the Ensembl ID, which could be annotated into the corresponding gene type using annotation information obtained from the Ensembl genome database (http://www.ensembl.org/). Ultimately, 11 types of lncRNAs were filtered from the multigene types. Noncoding RNA, lincRNA, processed transcript, TEC, bidirectional lncRNA, macro lncRNA, sense overlapping, sense intronic, retained intron, antisense lincRNA, and 3 prime overlapping ncRNAs were finally retained to generate a re-annotated lncRNA profile containing 911 lncRNAs.

Bioinformatics analysis

The GSE32894 series, which had a larger sample size (N=224), was used to determine the candidate prognostic lncRNAs. Through univariate Cox regression analysis of the GSE32894 series, lncRNAs that were significantly associated with disease-specific survival (DSS) were screened out as seed lncRNAs for further analysis (P<0.01). The random survival forest-variable hunting (RSF-VH) method was employed to identify the optimal prognosis-related lncRNAs. Here, we used the method with the value argument “nsplit” set to 10, “nrep” set to 100, and “nstep” set to 5, with 1000 trees grow and the k value set to 5 [18, 52]. The expression levels of the most valuable lncRNAs were compared among BLCA tissue (n=165), cancer-adjacent tissue (n=58) and normal bladder mucosa(n=9) in GSE13507.

Weighted by the regression coefficient, a risk score formula was constructed based on the expression level of each candidate lncRNA. Applying this formula, a risk score was calculated for every patient. By checking whether the risk score was greater than the median, the patients were separated into a high-risk (higher score) subgroup or a low-risk (lower score) subgroup. Kaplan-Meier estimates and the log-rank test were employed to assess the difference in prognosis between the two subgroups. Multivariate Cox regression analysis was conducted to determine predictive factors for BLCA prognosis and their independence from other clinical predictors. Using the “rms” package in R, two prognostic nomograms were established to predict the 3-year and 5-year DSS in GSE32894. One was based on the combination of the lncRNA signature with clinical factors, and the other was based on this signature alone. The predictive abilities of these nomograms were assessed with the concordance index (C-index) and calibration curves to compare the model-predicted values and actual observations of DSS. Receiver operating characteristic (ROC) curves were used to determine the prognostic value of the risk score, and the area under the curve (AUC) was also calculated for comparison with other clinical predictors using DeLong's test. In this method, patients having shorter DSS than the median DSS was labelled as positive and patients having longer DSS than the median DSS was labelled as negative. Patients surviving shorter than the median DSS time at the end of follow-up were excluded unless death had been observed. Correlation analysis using Pearson correlation coefficients was used to explore the associations among the 10 lncRNA-based risk scores and the expression levels of immune checkpoints and several potential therapeutic targets (P<0.01).

P < 0.01 was considered statistically significant in univariate Cox regression analysis and Pearson correlation analysis. In the rest of the methods, P<0.05 was considered statistically significant.

Gene set enrichment analysis (GSEA)

GSEA was conducted between the high-risk subgroup and low-risk subgroup from the GSE32894 series to identify potential cancer-related pathways. GSEA software V3.0 was employed to perform this analysis. Canonical representations of the biological process set, termed “c2.cp.v6.1.entrez.gmt” (1030 gene sets), were used in GSEA. The enrichment results were visualized in Cytoscape software V3.2.1. Dysregulated pathways were mapped in terms of biological processes. Gene sets with false discovery rate (FDR) values <25% and P<0.05 were termed “enriched” after performing 1000 random sample permutations.

ImmuCellAI analysis

ImmuCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) is an emerging tool used to estimate the abundance of 24 immune cells and overall infiltration scores based on a gene expression data set. The infiltrating data of the high-risk group in GSE32894 (N=112) and the low-risk group in GSE32894 (N=112) were obtained from the ImmuCellAI website. The abundance of each type of immune cell was then tested to detect the differences between the classified prognostic risk groups using the Wilcoxon test, and the results were visualized by a box plot. P<0.05 was considered statistically significant.