Here, we firstly illustrated the overall design and procedures of the research ( Figure 1 ).

We collected the somatic mutation profiles, transcriptome profiles and clinical information of all ccRCC patients recorded in the database of The Cancer Genome Atlas Program (TCGA, https://portal.gdc.cancer.gov/) on 16th September, 2020. The research included the data of all 539 ccRCC tumor samples. The data of somatic mutations contains the cumulative number, mutation site, and mutation type in each sample. The data of transcriptome profile contains the expression level of detected RNAs in each sample. From the clinical information, we selected age, gender, follow-up time, survival status, tumor grade, pathological stage, TNM stage for analyses (information of lymphatic metastasis was not analyzed as it is missing in most cases and is not routinely assessed as an important criterion for prognosis predicting in ccRCC). As short follow-up period usually means inaccurate conclusions, we discarded the cases of follow-up period less than 3 months in cox regressions and survival analysis concerning follow-up time.

The somatic mutation profiles of 336 cases were available and all cases were included in this part of analysis. Firstly, we calculated the cumulative number of somatic mutations for each patient. Then, we divided the patients into genomically stable (GS) group and genomically unstable (GU) group after ranking all patients according to the cumulative number of somatic mutations. Patients with top 25% numbers of somatic mutations were defined as GU group, and with last 25% numbers were defined as GS group. Finally, we performed differentially expressed (DE) analysis on all lncRNAs between the two groups. The analysis was conducted using Wilcoxon test with limma package (27) under R programming environment. The cut-off criteria for statistically significant difference were decided as false discovery rate (FDR) adjusted p < 0.05 and |log fold change (FC)| > 1. The identified DE-lncRNAs were regarded as GID-lncRNAs.

Of all GID-lncRNAs, overall survival (OS) related genes usually perform critical functions in the onset and progression of tumors, as well as make decisive difference on prognosis. Thus, we performed time-dependent univariate Cox regression with survival R package to extract OS-related GID-lncRNAs. The cut-off criterion for statistically significant correlations with OS was decided as p < 0.05. After discarding cases with follow-up period less than 3 months, 491 cases were included in this part of analysis.

The same 491 cases were included for constructing the risk signature. We divided all tumor cases into a training dataset and a testing dataset for developing and validating the risk signature. The cases were divided randomly according to a proportion of 7:3 (70% for training dataset and 30% for testing dataset, 377 cases in training dataset and 114 cases in testing dataset). We performed chi-square test to make sure that the training dataset and testing dataset were divided without bias.

Least absolute shrinkage and selection operator (LASSO) Cox regression was applied for develop a GID-lncRNAs based risk signature. The LASSO Cox regression would pick out variables for constructing the signature and give them coefficients. The GID-lncRNAs based risk signature was constructed as follow:

In both of the training dataset and the testing dataset, we computed the risk score for each patient according to the risk signature and then divided the patients into high-risk group and low-risk group from the median value of the risk score.

In order to assess the degree of GI and predict the prognostic outcomes of ccRCC patients with GID-lncRNAs, we have successfully developed the risk signature and divided patients into high-risk and low-risk groups. It is expected that patients in high-risk group are at GU condition. In order to confirm the supposition, we compared the mutation profiles of the two groups with maftools R package (28). All 336 cases with mutation profiles were included in this part.

After the development of the risk signature, we have to assess its reliability and robustness in estimating the prognostic outcomes of patients. The effectiveness of the signature was validated from the following 4 aspects:

LncRNAs do not code proteins themselves, but work through regulating protein-coding genes. In order to explore the potential function of the lncRNAs involved in the risk signature in ccRCC, we calculated the Pearson’s correlation coefficients between the lncRNAs and all mRNAs in the transcriptome profile to evaluate the co-expression relationships between lncRNAs and mRNAs. Then we extracted the top-300 related mRNAs for each lncRNAs and conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) (29) enrichment analysis on these related mRNAs, respectively. The cut-off criterion for significantly enriched items was decided as p < 0.05. The analysis was conducted using clusterProfiler R package (30).

In addition of predicting the prognostic outcome of ccRCC patients, the risk signature would be more serviceable if it could provide guidance for deciding treatment selections. With the development of precision medicine, ccRCC patients are given more opportunities and choices for curing or improving outcomes. We selected several targets of immunotherapy and targeted therapy approved by Food and Drug Administration (FDA) for analysis. The included targets of precise treatment are listed as follow: mTOR, KIT, PD-1, PD-L1, PDGFRA, PDGFRB, VEGFR1, VEGFR3, FLT3, RET, MET. We applied differentially expressed analysis for these targets between high-risk group and low-risk group to explore whether the expression level of these targets for precise treatments were related to different risk groups. The criteria for DE-genes were set as FDR adjusted p < 0.05 and |logFC| > 1.

Immunotherapy is a breakthrough for the treatment of ccRCC, but the overall response rate of PD-1/PD-L1 inhibitors is not satisfactory. It is certified that tumor mutation burden (TMB) could effectively predict the efficacy of immunotherapy of tumors as a biomarker (31). We calculated the TMB value (mutations per million bases) for each patient and compared the TMB values between high-risk and low-risk groups by t-test to explore whether the risk signature could help predict the response to immunotherapies targeting PD-1/PD-L1.

After computing the number of somatic mutations, patients were divided into GS-group and GU-group (The two groups contained 85 patients respectively, the ID and mutation number of patients are available in Table S1 ). We then performed DE analysis on all lncRNAs between GS-group and GU-group. We identified 46 DE-lncRNAs in all, which were regarded as functioning in GI (that is, GID-lncRNAs). Sixteen GID-lncRNAs were up-regulated in GU-group and 30 GID-lncRNAs were down-regulated in GU-group. The expression profiles of the GID-lncRNAs were shown in Figure 2A and the details (expression level, logFC, p value, FDR value) were given in Table S2 .

Time-dependent univariate Cox regression recognized 9 OS-related lncRNAs from all 46 GID-lncRNAs ( Figure 2B ). The official gene names of the 9 OS-related GID-lncRNAs were LINC00460, AL139351.1, AC156455.1, AL035446.1, LINC02471, LINC01234, AC022509.2, MIR222HG, and LINC01606.

Before developing the risk signature, we have divided all samples (with transcriptome profiles, follow-up time, and survival status) into training dataset and testing dataset. Chi-square test confirmed that there is no significant difference in clinical features between the two datasets. The information of patients in the two datasets and results of chi-square test are shown in Table 1 .

We developed the risk signature in the training dataset. All 9 OS-related lncRNAs were involved in LASSO Cox regression to construct the risk signature for predicting the prognostic outcomes of ccRCC patients. While constructing the risk signature, LINC01234 and MIR222HG were eliminated by LASSO Cox regression.

The GID-lncRNA based risk signature was given as follow:

The detailed information of the 7 lncRNAs was provided in Table 2 . Of the 7 lncRNAs composing the risk signature, 6 have positive coefficients and one has a negative coefficient. The up-regulation of lncRNAs with positive coefficients means worse outcome, that is, these lncRNAs may function as risky factors in ccRCC patients. Oppositely, lncRNA with negative coefficient may act as a protective factor.

For the subsequent analyses, we calculated the risk score for each patient with our risk signature and ranked patients according to their risk scores in training dataset and testing dataset, respectively. Then, patients were classified into high-risk group and low-risk group according to the median value of the risk score as a threshold (0.321 in training dataset and 0.329 in testing dataset). Samples in training dataset and testing dataset, their expression profile of the seven lncRNAs, as well as their risk scores and risk groups were given in Table S3 .

Figures 3A–C showed the landscape of the mutation profiles of ccRCC patients. According to different classified categories, missense mutation is the most frequently occurred type of mutations in ccRCC, and the amount of single nucleotide polymorphism (SNP) is significantly larger than that of insertion (INS) or deletion (DEL). VHL is the most frequently mutated gene in ccRCC, occupying almost 50% of all patients.

Figure 3D showed the mutation profiles of high-risk patients and low-risk patients separately. Mutations were more frequently occurred in high-risk group. Meanwhile, the mutation of several genes is rarely observed in low-risk group but frequently observed in high-risk group, for example, KDM5C. More types of mutations were observed in high-risk group. The analysis intuitively reflected the differences of mutation profiles between high-risk and low-risk groups.

Figure 5 exhibits the top-10 enriched KEGG-terms of KEGG enrichment on the co-expressed genes of the seven lncRNAs in the risk signature (Details and all terms of the results of KEGG enrichment analysis are available in Table S4 ). LINC00460 is closely related to signal pathways of “Cell cycle”, “DNA replication”, “p53 signaling pathway”, and “Mismatch repair”. AL139351.1 is involved in signal pathways concerning material metabolism and energy cycle. AC156455.1 is related to “RNA degradation”, “Drug metabolism”, and “Mismatch repair”. AL035446.1 takes part in signal pathways of “Cell cycel”, “p53 signaling pathway”, “DNA replication”, and “Drug metabolism”. LINC02471 mainly functions in signal pathways about material metabolism. AC022509.2 is related to “NF-kappa B signaling pathway”, “TNF signaling pathway”, “Apoptosis”, and “PD-L1 expression and PD-1 checkpoint pathway in cancer”. LINC01606 plays roles in “mTOR signaling pathway”, “Platinum drug resistance”, and “Drug metabolism”.

According to the results of KEGG enrichment analysis, we could conclude that the seven lncRNAs are involved in several biological themes closely related to the development and progression of tumors, such as cell cycle, DNA replication, and mismatch repair. Meanwhile, these lncRNAs may play important roles in material metabolism. Moreover, the functions of the lncRNAs have enriched in several widely admitted tumor-related signal pathways, such as NF-kappa B signaling pathway, p53 signaling pathway, mTOR signaling pathway, and TNF signaling pathway. Most importantly, the aberrant regulation of these lncRNAs may contribute to the resistance towards chemotherapy and immunotherapy.

DE analysis found the significantly different expression level of KIT and PD-1 between two risk groups ( Figure 6A ). PD-1 was up-regulated in the high-risk group (logFC = 0.55, FDR = 2.17E-5) and KIT was down-regulated in the high-risk group (logFC = -1.92, FDR = 1.34E-5). The results hint that patients in high risk group might response better to treatments targeting PD-1 and worse to treatments targeting KIT.

We calculated the TMB value for each patient and compared the inter-group difference of TMB value with t-test. The average TMB value in high-risk group was significantly higher than that in low-risk group. The TMB value was 1.25 in high-risk group and 1.00 in low-risk group (p = 0.019). Figure 6B exhibited the difference of TMB between two groups (one sample in high-risk group was omitted in the figure as its TMB value was too large). The results indicated a potential higher response rate towards immunotherapy in patients of the high-risk group classified by our risk signature.