A comprehensive prognostic and immunological analysis of telomere-related lncRNAs in kidney renal clear cell carcinoma

The expression of telomere-related genes and identification of TRLs

The study flow diagram is shown in (Supplementary Figure 1). We obtained the expression profiles of 2088 telomere-related genes telomere from 72 normal kidneys and 539 KIRC tissues. In KIRC, 37 telomere-related genes exhibited significant expression differences, including 15 down-regulated and 22 up-regulated genes (Figure 1A, 1B). The expression profiles of 37 differentially expressed telomere-related genes and lncRNAs were then utilized to conduct a Pearson correlation analysis. Under the conditions of |R| > 0.7 and P < 0.001, 13 telomere-related genes were retained and 648 lncRNAs were identified as TRLs. The co-expression networks between telomere-related genes and TRLs were illustrated in the Sankey diagram (Figure 1C).

Construction of TRLs prognostic signature

The univariate Cox analysis revealed that 15 of 648 TLRs were significantly related to overall survival (OS) (Figure 2A). In the training cohort, the aforementioned genes were subjected to the least absolute shrinkage and selection operator (LASSO) regression analysis for dimensionality reduction. There was a total of nine TRLs that were filtered out and used to create a prognostic signature (Figure 2B, 2C). The coefficient distributions of nine TRLs are illustrated in Figure 2D. The risk score for each KIRC patient can be estimated: risk score = (−0.387 × value of SIAH2-AS1) + (−0.218 × value of AL109741.1) + (−0.269 × value of LINC00551) + (0.097 × value of AC112715.1) + (0.379 × value of AC080013.2) + (−0.370 × value of LINC00571) + (−0.488 × value of AL161457.2) + (0.096 × value of AC002398.1) + (1.627 × value of AL512770.1). Additionally, we conducted a co-expression analysis of TRLs and telomere-related genes, and the results indicated that AC002398.1 and AL512770.1 co-expressed more than any other TRLs and telomere-related genes, which were linked to a poor prognosis (Figure 2E).

Validation of TRLs prognostic signature

KIRC patients were separated into low- and high-risk groups using the median score. Subsequently, we ranked the KIRC patients based on their risk scores from the training, testing, and entire cohorts (Figure 3A–3C). We then compared the patient survival status and the expressions of nine TRLs between two groups in the training, testing, and entire groups, respectively. We found that as the risk score increases, the number of deaths among patients with KIRC increased (Figure 3D–3F). Furthermore, four deleterious lncRNAs (AC112715.1, AC080013.2, AC002398.1, and AL512770.1) exhibited relatively high expression in the high-risk group, whereas five protective lncRNAs (SIAH2-AS1, AL109741.1, LINC00551, LINC00571, and AL161457.2) demonstrated the opposite trend (Figure 3G–3I).

The survival analysis revealed that the high-risk group suffered a significantly shorter OS rate than the low-risk group (Figure 4A–4C). In addition, there was a significant difference in PFS between the two groups, with the progression-free survival (PFS) rate being higher in the low-risk group (Figure 4D–4F). Except for stage N1, significant statistical differences in survival analysis were observed between KIRC patients with distinct clinical features (Supplementary Figure 2). The absence of a significant difference in the OS rate may be partially explained by the fact that there are fewer cases of N1 group. In addition, we found that patients with high-grade and advanced-stage tumors had significantly higher risk scorers than those with low-grade and early stages tumors (Supplementary Figure 3). Between these two groups, there was no significant difference in the distribution of age and gender. The preceding result revealed that the constructed prognostic signature holds for almost all clinical features of KIRC patients.

The predictive efficacy of the signature was further evaluated using the area under the curve (AUC) and its values for 1, 3, and 5 years were 0.673, 0.743, and 0.893, respectively (Figure 5A). Comparatively, the maximum AUCs for the testing cohort and the entire cohort were 0.686 and 0.753, respectively (Figure 5B, 5C). In addition, we observed that the risk score predicted the prognosis of KIRC patients better than other clinical variables (Figure 5D–5F).

Principal component analyses (PCA) were utilized to examine the distinct distribution patterns between two groups using a variety of data sets, including whole genome data, differentially expressed telomere-related genes, TRLs, and signature-based TRLs (Supplementary Figure 4). Using the expression of signature-based TRLs, the two groups tended to separate in two directions, whereas the other three expression profiles failed to distinguish between the groups. These results demonstrated that this signature possessed strong predictive ability.

The independence of TRLs prognostic signature

We utilized univariate and multivariate Cox analyses of variables to determine if the signature-based risk score was an independent prognostic factor for OS, independent of other clinical characteristics. In the training cohort, univariate Cox analysis revealed a significant association between risk score and patient survival (Figure 6A). In addition, a multivariate Cox analysis implied that the risk score was an independent prognostic factor (Figure 6B). The results of independent prognostic analysis for the testing (Figure 6C, 6D) the entire cohorts (Figure 6E, 6F) were identical. Accordingly, these findings suggested that the signature-based risk score was an independent prognostic factor for OS, suggesting which had a certain clinical utility for the prognosis evaluation. Subsequently, we developed a nomogram by integrating three independent indicators (age, grade, and stage) and a risk score to predict the 1, 3, and 5-years OS of KIRC patients (Figure 6G). The calibration curves validated the consistency of the nomogram (Figure 6H).

Functional enrichment analysis

Considering that the outcomes of KIRC patients in low- and high-risk groups differed significantly, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analyses were conducted to investigate the potential biological function of 534 genes with differential expression between the two groups. GO annotation revealed that these genes were significantly associated with biological behaviors like humoral immune response (Supplementary Figure 5A–5C). KEGG revealed that these genes were related to pathways including glycerophospholipid metabolism. (Supplementary Figure 5D, 5E). These biological behaviors and pathways are probably responsible for the high-risk group’s tendency to have poorer clinical outcomes.

Immune landscape of KIRC patients

Since the preceding results indicated a close association between risk scores and immune functions, we investigated the immune profiles of KIRC patients in low- and high-risk groups. A significant correlation between various immunocytes and the risk score was observed, which was determined by employing seven algorithms to study the relationship between risk score and tumor immune cell infiltration (Figure 7A). Then, using the single sample Gene Set Enrichment Analysis (ssGSEA) algorithm, we quantified scores for immune cell enrichment and function enrichment (Figure 7B). The proportions of CD8+ T cells, T follicular helper cell and tumor Infiltrating lymphocytes were higher in the high-risk group than the low-risk group (Figure 7C). As for the immunization-associated function enrichment scores, cytolytic activity, inflammation promoting, T cell co-inhibition/stimulation, and IFN response (Type 1) were higher in the high-risk group than the low-risk group (Figure 7D). We then determined that the immune score of the high-risk group was obviously higher than the low-risk group using the ESTIMATE algorithm (Figure 7E). Neither the stromal score nor the estimated score varied significantly between subgroups (Figure 7F, 7G). In addition, the poor prognosis was more prevalent among KIRC patients with higher immune scores (Figure 7H), indicating that the immune score can be used to predict the outcomes of these patients. However, there was no significant survival difference in the stromal score and estimate score (Figure 7I, 7J).

Tumor mutation analysis

To characterize tumor mutation analysis (TMB) in KIRC patients, we determined the TMB value for each individual case. The probability of a genetic mutation occurring in the low-risk group and high-risk group were 76.50% and 82.95%, respectively (Figure 8A, 8B). However, there is no significant difference in the value of TMB (Figure 8C). To investigate the value of TMB in KIRC, we then conducted survival analyses for different TMB groups and risk groups. The patient with a low TMB represented favorable OS than the patient with a high TMB (Figure 8D). The synergistic effect of TMB concentrations and risk scores on outcome prediction were then investigated in greater detail. Combining TMB subgroup analysis with the risk score revealed statistically significant differences in survivability across groups (Figure 8E). Additionally, the low-TMB combining low-risk group had the highest OS rate, whereas the high-TMB combining high-risk group had the lowest OS rate. These findings indicated that the signature-based risk score combined with TMB had a good prognostic value. Comparing the high-risk group to the low-risk group, the high-risk group exhibited a relatively higher immune status.

Immunotherapy responsiveness and drug sensitivity analysis

Subsequently, we determined the immune checkpoint inhibitors (ICIs) expression levels and drug sensitivity of KIRC patients. In the high-risk group, a number of ICIs, including PD-1, were upregulated (Figure 9A). Examining the differences in immunotherapy efficacy between the two groups, we then obtained the IPS of the KIRC patients from the Cancer Immunome Atlas (TCIA) database. Overall, there was no statistical difference in the IPS of CTLA4_negPD1_neg (Figure 9B); however, there existed difference in the group of CTLA4_posPD1_pos (Figure 9C), indicating that immunotherapy may be more suitable for the high-risk group. By comparing differences in the maximum 50% inhibition concentration (IC50) values of several drugs between two groups, we discovered that our signature may influence drug sensitivity to some extent. The fact that the IC50 values of 12 drugs in the low-risk groups were higher than those in the high-risk groups suggests that KIRC patients with high-risk scores may be more sensitive to the aforementioned chemical therapy (Figure 9D). The findings demonstrated that this prognostic signature has a certain guiding significance for selecting specific therapy for KIRC patients.

Data collection

The RNA sequencing data (FPKM value) of 539 KIRC and 72 normal patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Individuals with incomplete RNA expression data, missing prognostic information, or a 30-day survival time were excluded from our study. The 513 KIRC cases were randomly separated into two groups (1:1 ratio) for training (n = 257) and testing (n = 256). Patient demographics and baseline characteristics were listed in Supplementary Table 1.

Identification of telomere-related lncRNAs

TelNet is a user-friendly gene database designed to facilitate the study of genes involved in telomere maintenance mechanisms. Accordingly, the information on 2088 telomere-related genes were identified from TelNet (http://www.cancertelsys.org/telnet/) database [38]. With a preset threshold (|log2 fold change (FC)| ≥ 3.0 and false discovery rate (FDR) < 0.05), the “limma” package was used to identify the differentially expressed telomere-related genes [39]. The “limma” package was further used to conduct a Pearson correlation analysis on differentially expressed telomere-related genes and lncRNAs expression. When |correlation coefficient (R)| > 0.7 and P < 0.001, TRLs were eliminated from further consideration. A Sankey diagram was created to illustrate the relationship more intuitively between differentially expressed telomere-related genes and TRLs.

Construction and validation of the TRLs-based prognostic signature

We merged TRLs expression data with corresponding survival data and then conducted a univariate Cox analysis to identify TRLs associated with OS (P < 0.05). In the training cohort, the LASSO regression analysis was used to develop the optimal prognostic signature. The risk score was estimated in the three cohorts using the following equation: $\text{Risk score}={\displaystyle {\sum }_{i=1}^{n}x(i)\times Coef(i).}$ The x(i) and Coef(i) represented the expression values of TRLs and their coefficients, respectively. Using the median score, the patients were classified into low-risk (< median) and high-risk (≥median) groups. Additionally, Kaplan-Meier survival analysis was utilized to compare the OS and PFD between the two groups using log-rank tests. Moreover, the receiver operating characteristic (ROC) curves and subgroup analyses were conducted to further evaluate the signature’s ability to predict outcomes and the performance of relevant clinical parameters. To determine the independence of signature-based risk scores, Cox regression analyses were conducted. The aforementioned analyses were validated in the testing and entire cohorts. A nomogram was developed eventually by integrating all independent prognostic parameters to qualitatively predict the OS of KIRC patients across the entire cohort.

Functional enrichment analysis

To investigate the potential molecular mechanisms underlying the prognostic signature, the differentially expressed genes between low- and high-risk groups were identified (|log2FC| > 1 and FDR < 0.05) using the “limma” package. The functional enrichment analysis of GO and KEGG was then carried out on the above genes. An adjusted P < 0.05 was regarded to be the cut-off point.

Analysis of tumor immune microenvironment

Seven algorithms (CIBERSORT-ABS, CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER, XCELL) were used for calculating the immune cell infiltration in various KIRC cases. The correlation between infiltration levels of immune cell and risk score was then determined using Spearman correlation analysis. In addition, the ssGSEA was employed for analyzing the relative abundance of immune populations between different groups [40]. Subsequently, we used the ESTIMATE algorithm to obtain the immune score for each KIRC patient and for comparing the score differences between risk groups [41].

Tumor mutation analysis

Based on somatic mutation frequencies, TMB was calculated. The median TMB value was then used to divide the KIRC patients into low- and high-TMB groups. The waterfall plot of a mutational landscape was visualized through the “maftool” package [42]. In addition, comparative and survival analyses were conducted between two groups to investigate the difference in somatic mutation.

Prediction of immunotherapy responsiveness and potential drug sensitivity

We also compared the expression levels of several ICIs genes between low- and high-risk groups and then utilized the TCIA database to predict potential immunotherapy responses [43]. Gene expression profiles and drug response data were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database. Additionally, we adopted ridge regression to predict the IC50 of drugs with the “oncoPredict” package [44]. The IC50 was calculated as tolerance ability, and a higher IC50 indicated the cells showing higher resistance to drugs. The difference in sensitivity score was then analyzed between the two groups.

Statistical analysis

All statistical analyses were performed in R statistical software (version 4.0.1) and Perl language (version 5.30.2). Unless otherwise specified, statistical significance was considered for two-tailed P < 0.05.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.