The FPKM-RNA sequence and clinical information of ccRCC were downloaded from the TCGA-KIRC data portal (https://portal.gdc.cancer.gov/), where it contained a total of 537 ccRCC tissues and 72 adjacent tissues. The clinical information mainly included age, gender, clinical stage, T stage, N stage, M stage, survival status, survival time, and survival prognosis. The exclusion criteria were as follows: 1) The pathological diagnosis did not meet ccRCC; 2) The RNA sequence and clinical data were incomplete; and 3) The follow-up time did not exceed 30 days. Subsequently, the expression data of these RNAs were sorted, annotated, and then assigned to protein-coding genes and lncRNAs using the Ensembl human genome browser (http://asia.ensembl.org/info/data/index.html) and Perl program. The extracted data were normalized and processed by log2 transformation. The expression data of 89 ccRCC patients from the ICGC database (https://dcc.icgc.org/analysis) were obtained for the external validation of the FRlncRNA signature. The “limma” package in R software was utilized to correct the transcriptome data we have downloaded.

The FerrDb database (http://www.zhounan.org/ferrdb/legacy/index.html) was used to obtain the FRG dataset, containing a total of 214 FRGs (Supplementary data) where 203 FRGs were found in the TCGA dataset. Among them, 62 FRGs were differentially expressed in ccRCC. The relation between FRGs and lncRNAs was analyzed by the Pearson correlation analysis. LncRNAs would be included in this study when the square of the correlation coefficient (R ²) was greater than 0.3, and concomitantly, the p-value was lower than 0.001. Finally, the “limma” package in R software was applied to extract a total of 1,669 FRlncRNAs (Supplementary data).

To establish an effective prognostic prediction model, we randomly divided 537 ccRCC patients into training cohorts and testing cohorts in a 1:1 ratio. Finally, 243 people were included in the training cohorts, and 264 people were included in the testing cohorts, according to the exclusion criteria. The detailed patient grouping process is shown in Supplementary Figure S1. The basic characteristics of each group are showed in Table 1 (Details in Supplementary data). The FRlncRNA signature was built according to the training cohorts, and the capability of predicting prognosis was evaluated via the testing cohorts, overall cohorts, and ICGC cohorts. The prognostic capability of FRlncRNAs in the training cohorts was assessed by the univariate Cox regression analysis. If p < 0.05, it would be included in the least absolute shrinkage and selection operator (Lasso) regression using the “glmnet” package in R software to avoid overfitting. The risk score of each patient, established by incorporating the Lasso regression into the multivariate Cox regression analysis, was calculated according to the function, ∑ⁿ _i=1β_i∗ (expression of lncRNA_i), where β represented the regression coefficient. Patients were classed into high- and low-risk groups based on the median risk score, and the survival rate between two groups was compared using the log-rank test.

The risk score of each patient in the test cohorts, overall cohorts, and ICGC cohorts was calculated using the same way to confirm the stability of the established model. The survival outcomes of each cohort were analyzed by the Kaplan–Meier (K-M) survival curve. The “ROC package” in R software was employed to analyze the specificity and sensitivity of the established model based on the ROC curve and its area under the curve (AUC) value.

A prognostic nomogram according to the aforementioned risk score and traditional prognosis-related clinical variables (age, grade, and stage) was established to statistically predict the prognosis of ccRCC patients. Subsequently, the reliability and accuracy of the nomogram were evaluated via the concordance index (C-index), calibration curve, and ROC curve. The basic characteristics of patients were included into the multivariate Cox regression analysis to determine whether the risk score was an independent predictor of prognosis.

The gene set enrichment analysis (GSEA) software was employed to perform the explanation of the functional enrichment of these FRlncRNAs. The pathway components of the high-risk group and the differences of the pathway activity and expression patterns, which were downloaded from MSigDB and GSEA 4.1.0, were analyzed in the pathway analysis dataset including c2. cp.kegg. v7.4. symbols and c5. go.bp. v7.4. symbols. A two-tailed p-value less than 0.05 was considered to be significant. To clarify how target genes of these FRlncRNAs participated in the development of ccRCC, we performed a functional analysis of related FRGs. The “clusterProfiler” package (Yu et al., 2012) and “org.Hs.eg.db” package in R software were used for functional enrichment of target genes based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway and gene ontology (GO) enrichment analysis. The cut-off criterion was set as p less than 0.05 and a q value more than 0.05.

When the co-expression coefficient was greater than 0.3 and the p-value was less than 0.001 using the “limma” package, we believed that there was a good correlation between FRlncRNAs and FRGs. Cytoscape 3.6.0 software was used to visualize the network between eight FRlncRNAs and related FRGs. After that, we further identified the major FRlncRNAs and related FRGs by searching the relevant literature and multiple databases. We also explored the ceRNA network of major FRlncRNAs. The gene expression profiling interactive analysis (GEPIA, http://gepia.cancer-pku.cn/) contained RNA-seq and clinical data compiled by TCGA and GTEx after standardized analysis. The UALCAN online database (http://ualcan.path.uab.edu/index.html) reported the differences of the target gene expression between normal tissues and ccRCC-graded tissues. The K-M Plotter database (http://kmplot.com/analysis/) included data on the correlation between the gene expression and prognostic data of 530 patients with ccRCC. The expression levels of the major FRlncRNAs-FRGs and the prognostic correlation were verified in the above three databases.

A total RNA extraction micro kit (RNT411-03, Mabio, Guangdong, China; http://www.mabiotech.cn/en/list-643-1.html) was used to extract the total RNA from tumor tissues and normal tissues based on the operation instructions, and then, a spectrophotometer was employed to detect the concentration and check the quality of total RNA. Then, a cDNA synthesis kit with random primers (AG11711, AG, Changsha, China; https://agbio.com.cn/product/evo-m-mlv-rt-kit-with-gdna-clean-for-qpcr-ii/?v=b838b393d55f) was used in a 20 μ1 reaction volume with 1 μg total RNA for cDNA synthesis. The mRNA expression was detected by the Script SYBR Green PCR kit (AG11702, AG, Changsha, China; https://agbio.com.cn/product/sybr-green-premix-pro-taq-hs-qpcr-kit-ii/?v=b838b393d55f) via an ABI7900HT Fast Real-time PCR machine, with the following conditions: pre-denaturation at 95°C for 30 s, followed by 40 cycles of denaturation at 95°C for 5 s, and annealing and extension at 60°C for 30 s. GAPDH was the internal control for mRNA, and the 2-ΔΔCt function was used to calculate the gene expression level (DiMagno, 2007). QRT-PCR primer sequences used in our experiments are listed as follows: GAPDH forward primer (F): 5′-TGA​CTT​CAA​CAG​CGA​CAC​CCA-3’; GAPDH reverse primer (R): 5′-CAC​CCT​GTT​GCT​GTA​GCC​AAA-3’; LINC00460 F: 5′-TAA​ACC​TAG​GGG​CCG​GTC​G-3’; LINC00460 R: 5′-AAC​GGT​CCA​GAG​CAG​ACA​AAA-3’; LINC00944 F: 5′-AGA​CGC​ACA​TCA​GGA​AGA​CAG-3’; LINC00944 R: 5′-TGA​GTT​ACA​GGG​ACC​GAA​GC-3’; LINC01550 F: 5′-GGT​GCA​GTC​TCC​TCA​GAA​CTA​C-3’; LINC01550 R: 5′-GGG​AGA​GGG​AGA​ACG​ACT​GT-3’; EPB41L4A-DT F: 5′-CGG​AGC​AGG​TGC​AAT​CTG​T-3’; and EPB41L4A-DT R: 5′-TCA​AAA​CTA​CGT​CTG​ATG​CCA​AA-3’.

The two-tailed Student’s t-test was calculated for the difference within groups or among groups. Categorical variables were presented as proportions, and the chi-square test was used for the comparisons between groups. Multivariate or univariate Cox proportional hazard regression analysis was calculated for evaluating the prognostic significance. The log-rank test and K-M curve assessed the prognostic results. R software (version 3.6.0) was used to draw the heatmap, GSEA, survivorship curve, ROC curve, nomogram, and calibration plot. A two-tailed p-value less than or equal to 0.05 was considered as the statistical significance.

The flowchart of this work is showed in Figure 1. The univariate Cox regression analysis was used to analyze the expression of FRlncRNAs in the training cohorts. A total of 678 lncRNAs were found to be closely related with the prognosis of ccRCC. High overfitting of these prognostic-related lncRNAs were eliminated by Lasso Cox analysis, and 15 FRlncRNAs were identified (Figure 2A; Supplementary data). Then, multivariate Cox regression analysis extracted a prognostic signature containing eight FRlncRNAs and their coefficients (Supplementary Table S1; Figure 2B), calculated with the following formula: Risk score = (0.190*AL590094.1) + (0.047*LINC00460) + (0.204* LINC00944) + (0.168*AC024060.1) + (0.116*HOXB-AS4) + (0.048* LINC01615)—(0.134* EPB41L4A-DT)—(0.345* LINC01550). Based on the hazard ratio (HR) score gained by the multivariate Cox regression analysis, AC024060.1, AL590094.1, HOXB-AS4, LINC00460, LINC00944, and LINC01615, whose HRs were greater than 1, were risk factors, but EPB41L4A-DT and LINC01550, whose HRs were less than 1, were protective factors (Supplementary Table S1).

To evaluate the sensitivity and stability of the prognostic risk score, the training cohorts were divided into the low-risk group (122 cases) and the high-risk group (121 cases) based on the median of risk scores (0.86). The results of the K-M curve showed that the survival ability of patients in the high-risk group was significantly lower than that in the low-risk group (p < 0.001) (Figure 2C). The accuracy of the prognostic signature was evaluated by the ROC curve, and the AUC value was 0.875 (Figure 2D). The AUC values of 1, 3, 5, and 10 years of overall survival (OS) were 0.759, 0.818, 0.871, and 0.914, respectively (Figure 2E). The heatmap showed remarkable differences in the expression of eight FRlncRNAs between the high-risk group and the low-risk group (Figure 2F). The scatter plot indicated that ccRCC patients with a high risk score had a lower survival rate than those with a low-risk score (Figure 2G). Moreover, the distribution map of the risk score was consistent with the categorization of patient groups (Figure 2H). The prognostic effect of eight FRlncRNAs evaluated by K-M curves displayed that higher expressions of AC024060.1, AL590094.1, HOXB-AS4, LINC00460, LINC00944, and LINC01615 and lower expressions of EPB41L4A-DT and LINC01550 were linked to inferior OS (p < 0.01) (Figures 3A–H). The prognostic risk-related model we constructed exhibited a good stability and sensitivity in predicting the OS of ccRCC patients.

To validate the predictive capacity of the FRlncRNA signature, risk scores of patients were calculated in the testing cohorts and overall cohorts, and patients were classified into the low-risk group and the high-risk group based on the median of risk scores (0.85 and 0.86, respectively). The OS in the testing cohorts (p < 0.001) (Figure 4A) and overall cohorts (p < 0.001) (Figure 4C) were analyzed by K-M curves, demonstrating that these results were in line with the training cohorts. The ROC curves of testing cohorts (AUC = 0.690) (Figure 4B) and overall cohorts (AUC = 0.771) (Figure 4D) illustrated that the FRlncRNA signature has an accurate predictive capability about the OS of ccRCC patients, which was further validated by ROC time curves and their AUC values, such as the AUC value in testing cohorts, 1-year AUC = 0.757, 3-year AUC = 0.684, 5-year AUC = 0.681, and 10-year AUC = 0.725 (Figure 4E) and the AUC value in overall cohorts (1-year AUC = 0.758, 3-year AUC = 0.750, 5-year AUC = 0.764, and 10-year AUC = 0.812) (Figure 4G). The consistent expression profiles of eight FRlncRNAs in the training cohorts are shown in the heatmap (Figures 4F,H). A lower survival rate was observed in the high-risk group than the low-risk group, and distribution maps of the risk score validated a higher risk score in the high-risk group (Figures 4I–L). In addition, ICGC cohorts were used to evaluate the constructed model, which showed a good ability to predict the survival rate of patients with ccRCC (Supplementary Figure 2A–E). These results showed that the FRlncRNA signature can be a good indicator in predicting the prognosis of patients compared with other existing signatures reported in recent studies (Supplementary Table S2) (Canxuan and Dan, 2021; Hong et al., 2021; Ma et al., 2021; Xing et al., 2021; Yu et al., 2021; Zheng et al., 2021). Taken together, our data implied that the FRlncRNA signature showed a stable prognostic-predictive ability.

The stratified analysis of clinicopathological characteristics was performed to assess the predictive ability of the FRlncRNA signature and the stability of its OS prediction in the high-risk and low-risk groups including gender (female and male), age (≤65 years old and >65 years old), grade (I-II and III-IV), stage (I-II and III-IV stage), T stage (T1-2 stage and T3-4 stage), and M stage (M0 stage and M1 stage). The results of the K–M curve in different clinical characteristics suggested that the OS of the high-risk group was worse than that of the low-risk group (p < 0.01) (Figures 5A–L).

The risk score was demonstrated by the univariate and multivariate Cox regression analyses to be an independent prognostic factor (p < 0.05) (Supplementary Table S3, Figures 6A,B). Subsequently, clinicopathological characteristics including the age, grade and stage, and risk score were applied to construct the nomogram using the “rms” package in R software to predict the 1-, 3-, and 5-year OS of ccRCC patients (Figure 6C). The increasing risk score deteriorated the prognosis of ccRCC. The results of the multivariate ROC curve suggested that the AUC value was 0.771, which was higher than that of the grade (0.664) and stage (0.714), implying that the nomogram had the ability of accurate prediction for survival outcomes of ccRCC (Figure 6D). We used the calibration curve to observe whether the actual prognostic value was consistent with the predicted value of the nomogram and found that the calibration curves of 1-, 3-, and 5-year survival rates were consistent with the nomogram (Figure 6E). The clinical influences of the risk score for ccRCC patients in the training, testing, and overall cohorts are showed in Supplementary Table S4.

The underlying molecular mechanisms of the FRlncRNA signature involved in the high-risk group were further verified by GSEA . Signaling pathways including the P53 signaling pathway [enrichment score (ES) 0.49; normalized enrichment score (NES) 1.76; nominal (NOM) p-value 0.03], the cytokine–cytokine receptor interaction signaling pathway (ES 0.39; NES 1.71; NOM p-value 0.03), the tumor necrosis factor–mediated signaling pathway (ES 0.50; NES 1.87; NOM p-value 0.005), and regulation of the T helper 1 type immune response signaling pathway (ES 0.64; NES 1.98; NOM p-value 0.007) were markedly enriched in the high-risk group (Figure 7).

The top 30 terms from the GO analysis of FRlncRNA-related FRGs are demonstrated in the dot plot (Figures 8A–C). GO analysis consisted of biological process (BP) analysis mainly including positive regulation of the catabolic process and intrinsic apoptotic signaling pathway; cellular component (CC) analysis mainly containing focal adhesion, cell-substrate adherens junction, and cell-substrate junction; and molecular function (MF) analysis mostly composing of protein serine/threonine kinase activity, ubiquitin protein ligase binding, and iron ion binding. The “pathway-gene network” and “pathway-gene clustering,” as shown in Figures 8D–F, were plotted to represent the complex relationship between FRlncRNA-related FRGs and KEGG pathways.

The co-expression network between eight FRlncRNAs and 49 FRGs (R ² > 0.3 and p < 0.001) is shown in Figure 9A. The Sankey diagram showed the interrelation between eight FRlncRNAs, 49 FRGs, and the risk type (Figure 9B). By analyzing the results of differential expression analysis (p < 0.01) and searching relevant literatures (Shijie Zhang et al., 2021; Chen and Zheng, 2021; Xuan et al., 2021) and databases, we selected four FRlncRNAs (LINC00460, LINC00944, LINC01550, and EPB41L4A-DT) for further study. Compared with normal tissues, LINC00460 (logFC = 5.039, p = 1.74E-19) and LINC00944 (logFC = 3.906, p = 2.46E-32) were significantly upregulated in ccRCC tissues, and LINC01550 (logFC = 0.359, p = 0.005) and EPB41L4A-DT (logFC = 0.2400, p = 0.003) were significantly downregulated in ccRCC tissues. To more accurately predict their action pathways, the ceRNA network (lncRNA–miRNA–mRNA) of four FRlncRNAs was explored (Supplementary Figure S3).

Two databases (GEPIA and K-M Plotter) were employed to investigate four FRlncRNAs including expression levels and survival outcomes. The expression levels of LINC00460 (Pr = 7.18E-09) and LINC00944 (Pr = 2.54E-11) increased gradually with the increase of stages but those of LINC01550 (Pr = 0.000122) and EPB41L4A-DT (Pr = 1.06E-10) decreased gradually with the increase in stages, suggesting that four FRlncRNAs were strongly correlated with tumor progression (Figures 10A–H). High expression levels of LINC00460 and LINC00944 and low expression levels of LINC01550 and EPB41L4A-DT were related to the poor prognosis (p < 0.05) (Figures 10I–P). Similarly, lower 5-year OS (p < 0.05) was noted in 530 ccRCC patients from the K-M Plotter database with an increasing expression of LINC00460 or decreasing expressions of LINC01550 and EPB41L4A-DT (Figures 10Q–S).

Through linear correlation analysis (R > 0.3 and p < 0.001) and literature retrieval (Shao et al., 2019; Xiong et al., 2021a; Hong et al., 2021), we selected three target genes (BNIP3, RRM2, and GOT1) for further verification. BNIP3 (p < 0.05) and RRM2 (p < 0.05) were significantly upregulated in ccRCC tissues, but GOT1 (p < 0.05) was significantly downregulated in ccRCC tissues (Figures 11A–C). The expression levels of BNIP3 and GOT1 decreased gradually with the increase in grades, but RRM2 increased gradually with the increase in clinical grades and stages, which revealed that three target genes were closely related to tumor progression (Figures 11D–I). High expression levels of RRM2 and low expression levels of BNIP3 and GOT1 were related to poor prognosis after the analysis of GEPIA and K-M Plotter databases (p < 0.05) (Figures 11J–O). The linear correlation analysis found that RRM2, BNIP3, and GOT1 may be potential targets of LINC00460 (R = 0.35; p < 2.2E-16), LINC01550 (R = 0.57; p < 2.2E-16), and EPB41L4A-DT (R = 0.33; p < 3.8E-15), respectively, which also needed further experimental verification (Figures 11P–R).

The expression levels of four FRlncRNAs were verified by qRT-PCR in the tumor and adjacent normal tissues collected from twenty ccRCC patients (Supplementary Table S5). LINC01550 and EPB41L4A-DT in tumor tissues were downregulated, while LINC00460 and LINC00944 were upregulated, showing the statistical significance (Figures 12A–D, t-test, p < 0.05; Figures 12E–H, paired t-test, p < 0.05). These results were consistent with our previous verification results, but more samples were still needed for verification.

We summarized the result of 507 ccRCC patients calculated by various algorithms and compared all the immune cell subtypes in two groups. Infiltration proportion of partial cell subtypes had an obvious difference between two groups, among which mainly T cell NK, B cell, T cell follicular helper, and T cell regulatory (Tregs) had a higher infiltration proportion in the high-risk group, while T cell CD4⁺, neutrophils and endothelial cells had a lower proportion (Figure 13A). The immune functions of the high-risk group and low-risk group were analyzed, respectively, using “GSVA” and “GSEABase” packages in R software and found that almost all items (APC co-stimulation, CCR, check-point, cytolytic activity, inflammation promoting, para inflammation, T-cell co-inhibition, T-cell co-stimulation, Type I IFN response) in the high-risk group were upregulated (p < 0.05, Figure 13B), indicating a significant change in the immunophenotype in the high-risk group. We further explored the expression of immune checkpoint-related markers in two groups and found some markers (CTLA4, CD40LG, LAG3, CD44, CD27, CD160, TNFRSF18, CD40, TNFSF4, CD244, TMIGD2, LAIR1, TIGIT, TNFRSF9, PDCD1, CD86, IDO2, CD200R1, BTLA4, TNFSF9, LGALS9, CD70, CD48, CD80, TNFRSF25, ICOS, TNFSF14, CD28, and TNFRSF8) in the high-risk group were upregulated, and some markers (NRP1, KIR3DL1, HHLA2, TNFSF18, HAVCR2, and TNFSF15) were downregulated, indicating an immunosuppressive and exhausted phenotype in the high-risk group (Figure 13C). Based on the above analyses, we found two groups had a significant distinct pattern of immune infiltration, which may lead to different survival benefits.