Advertisement

Establishing a Prognostic Model Based on Three Genomic Instability-related LncRNAs for Clear Cell Renal Cell Cancer

  • Author Footnotes
    † These authors have contributed equally to this work
    Shen Lulu
    Footnotes
    † These authors have contributed equally to this work
    Affiliations
    Department of Obstetrics and Gynecology, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Author Footnotes
    † These authors have contributed equally to this work
    Hou Hualing
    Footnotes
    † These authors have contributed equally to this work
    Affiliations
    Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Zhang Shan
    Affiliations
    Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Chen Dianxi
    Affiliations
    Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Li Yiqing
    Correspondence
    Address for correspondence: Li Yiqing, Li Qin, MD, Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, 1277 Jiefang Avenue, Wuhan, Hubei Province, 430022, P.R. China
    Affiliations
    Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Li Qin
    Correspondence
    Address for correspondence: Li Yiqing, Li Qin, MD, Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, 1277 Jiefang Avenue, Wuhan, Hubei Province, 430022, P.R. China
    Affiliations
    Department of Vascular Surgery, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
    Search for articles by this author
  • Author Footnotes
    † These authors have contributed equally to this work
Published:February 23, 2022DOI:https://doi.org/10.1016/j.clgc.2022.02.005

      Abstract

      Background

      Among all types of renal cell cancer (RCC), clear cell renal cell cancer (ccRCC) is the most common and aggressive one. Emerging evidence uncovers that long non-coding RNAs (lncRNAs) are involved in genomic instability, which correlates to the clinical outcomes of patients who suffer from various kinds of cancers.

      Methods

      We gathered expression profiles of transcriptome RNA and clinical information about ccRCC from The Cancer Genome Atlas (TCGA) and The Gene Expression Omnibus (GEO) database. The lncRNA expression profiles and somatic mutation data were combined to identify genome instability-related lncRNAs (GILncRs) by significance analysis of T test. By means of univariate and multivariate cox regression analyses, 3 GILncRs strongly associated with patient prognosis were screened out to build a genomic instability-related risk score (GIRS) model. We use R-version 4.0.4 to draw Kaplan-Meier plots and ROC curves for survival prediction.

      Results

      The somatic mutation count was higher in genomic unstable group. PBRM1 showed lower expression in genomic unstable group. Three lncRNAs such as LINC00460, AC156455.1, LINC01606 were included in the GIRS model. Patients had poorer prognosis with higher risk score of GIRS model. The somatic mutation count was higher in patients with higher risk score while PBRM1 expression was lower. The GIRS model was independent from other clinical factors. The GIRS model was superior to other 2 published lncRNA signatures in survival prediction.

      Keywords

      Introduction

      As is well known that, the most common pathological subtype of kidney cancer is renal cell cancer (RCC), which contains various types of cancer. In 2004, The World Health Organization (WHO) identified 11 histologic types of RCC including clear cell, chromophobe, oncocytoma and papillary of 2 types, etc.
      • Capitanio U.
      • Bensalah K.
      • Bex A.
      • et al.
      Epidemiology of renal cell carcinoma.
      ,
      • O'Rourke C.J.
      • Knabben V.
      • Bolton E.
      • et al.
      Manipulating the epigenome for the treatment of urological malignancies.
      Clear cell renal cell cancer (ccRCC) is the most common and the most aggressive subtype of RCC.
      • Linehan W.M.
      Genetic basis of kidney cancer: role of genomics for the development of disease-based therapeutics.
      Moreover, one-third of ccRCC patients with locally advanced or distant spread are categorized into late-stage tumor patients, which results in extremely poor prognosis.
      • Vuong L.
      • Kotecha R.R.
      • Voss M.H.
      • Hakimi A.A.
      Tumor microenvironment dynamics in clear-cell renal cell carcinoma.
      Additionally, traditional treatments do not work since the ccRCC patients are not sensitive to both chemo- and radio-therapy,
      • Rydzanicz M.
      • Wrzesinski T.
      • Bluyssen H.A.
      • Wesoly J.
      Genomics and epigenomics of clear cell renal cell carcinoma: recent developments and potential applications.
      while immunotherapy such as immune checkpoint inhibitors (ICIs) have been confirmed to be highly efficacious, in which CTLA-4, PD-1/PD-L1 are examples.
      • Song M.
      Recent developments in small molecule therapies for renal cell carcinoma.
      • Sanchez-Gastaldo A.
      • Kempf E.
      • Gonzalez Del Alba A.
      • Duran I.
      Systemic treatment of renal cell cancer: A comprehensive review.
      • Guo C.
      • Zhao H.
      • Wang Y.
      • et al.
      Prognostic value of the neo-immunoscore in renal cell carcinoma.
      Genomic instability (GI) represents a notably important feature of cancer, which results from mutations in DNA repair genes
      • Negrini S.
      • Gorgoulis V.G.
      • Halazonetis T.D.
      Genomic instability–an evolving hallmark of cancer.
      and subsequently promotes tumor growth and metastasis as well as resists to various therapies.
      • Kalimutho M.
      • Nones K.
      • Srihari S.
      • Duijf P.H.G.
      • Waddell N.
      • Khanna K.K.
      Patterns of genomic instability in breast cancer.
      Accumulating evidence shows that genomic instability correlates to the clinical outcomes of patients suffered from various cancers.
      • Kalimutho M.
      • Nones K.
      • Srihari S.
      • Duijf P.H.G.
      • Waddell N.
      • Khanna K.K.
      Patterns of genomic instability in breast cancer.
      • Poole L.A.
      • Cortez D.
      Functions of SMARCAL1, ZRANB3, and HLTF in maintaining genome stability.
      • Syed A.
      • Tainer J.A.
      The MRE11-RAD50-NBS1 complex conducts the orchestration of damage signaling and outcomes to stress in DNA replication and repair.
      Long non-coding RNAs (lncRNAs) pertain to non-coding RNAs (ncRNAs) referred to as those RNA transcripts which are lack of protein-coding potential. If the cut-off value of the length > 200 nucleotides,
      • Cheetham S.W.
      • Gruhl F.
      • Mattick J.S.
      • Dinger M.E.
      Long noncoding RNAs and the genetics of cancer.
      ,
      • Moghaddas Sani H.
      • Hejazian M.
      • Hosseinian Khatibi S.M.
      • Ardalan M.
      • Zununi Vahed S.
      Long non-coding RNAs: An essential emerging field in kidney pathogenesis.
      the ncRNAs are defined as long ncRNAs (lncRNAs), otherwise they are detected to be small ncRNAs.
      • Martens-Uzunova E.S.
      • Bottcher R.
      • Croce C.M.
      • Jenster G.
      • Visakorpi T.
      • Calin G.A.
      Long noncoding RNA in prostate, bladder, and kidney cancer.
      Emerging evidence indicates that ncRNAs are crucial for DNA repair and genomic integrity, which plays an important role in survival of tumor cell and prevention of tumorigenesis. As reported by Mathias et al, NORAD, which belongs to highly conserved lncRNAs, is essential for maintaining genomic stability since it is an indispensable part of the topoisomerase complex NARC1. Moreover, NORAD is upregulated due to DNA damage, which induces chromosomal instability.
      • Lee S.
      • Kopp F.
      • Chang T.C.
      • et al.
      Noncoding RNA NORAD regulates genomic stability by sequestering PUMILIO proteins.
      ,
      • Munschauer M.
      • Nguyen C.T.
      • Sirokman K.
      • et al.
      The NORAD lncRNA assembles a topoisomerase complex critical for genome stability.
      As reported by Wu et al, knocking out the lncRNA GUARDIN, which is tightly related to p53 response, results in inhibiting tumor growth and increasing cytotoxicity induced by additional stress. Thus, GUARDIN is important for maintaining genomic integrity and may contribute to the development of targeted therapy.
      • Hu W.L.
      • Jin L.
      • Xu A.
      • et al.
      GUARDIN is a p53-responsive long non-coding RNA that is essential for genomic stability.
      However, clinical significance and prognostic value of lncRNAs related to genomic instability still remain unexplored in spite of the above findings.
      • Khanduja J.S.
      • Calvo I.A.
      • Joh R.I.
      • Hill I.T.
      • Motamedi M.
      Nuclear noncoding RNAs and genome stability.
      • Khorkova O.
      • Hsiao J.
      • Wahlestedt C.
      Basic biology and therapeutic implications of lncRNA.
      • Rothschild G.
      • Basu U.
      Lingering questions about enhancer RNA and enhancer transcription-coupled genomic instability.
      • Wan G.
      • Liu Y.
      • Han C.
      • Zhang X.
      • Lu X.
      Noncoding RNAs in DNA repair and genome integrity.
      Here, we attempted to figure out genomic-instability associated lncRNAs and establish a genomic instability-related risk score model to predict the clinical outcomes of ccRCC patients.

      Materials & Methods

      Data Collection and Pretreatment

      We downloaded the data of transcriptome RNA-sequencing expression, somatic mutation and clinical information for ccRCC patients from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), which were gathered from 539 tumor tissues. We also downloaded data from The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) with accession number GSE22541. The patients whose overall survival (OS)≤ 30 days were excluded as they possibly died for unforeseen reasons like cardiovascular accident and infection. With a use of the merge script in Perl language (http://www.perl.org/), the RNA-sequencing expression data were integrated into a matrix file and each RNA was identified as mRNA or lncRNA. The Ensembl database (http://asia.ensembl.org/index.html) translated the Ensembl IDs to gene symbols.

      Identification of Genomic Instability-Related lncRNAs (GILncRs)

      To identify the genomic instability-related lncRNAs (GILncRs) of ccRCC, lncRNA expression data and somatic mutation information were merged into a matrix file and then processed following the steps below: (1) The cumulative somatic mutation count of each patient was calculated; (2) all patients were sorted in increasing order of the cumulative somatic mutation count; (3) the top 25% patients were considered as low mutation (LM) group, and the bottom 25% were considered as high mutation (HM) group; (4) use the significance analysis of the T test to compare the expression levels of lncRNA between the LM group and the HM group; (5) differently expressed lncRNAs (false discovery rate adjusted P < .05) were defined as GILncRs. Depend on the expression level of GILncRs, all patients were divided into genomic unstable (GU) group with high somatic mutation or genomic stable (GS) group with low somatic mutation using the cluster analysis. The heat map was drawn to show the difference of GILncRs expression between GU group and GS group.

      Functional Enrichment Analysis

      To evaluate the relationship between the expression of lncRNAs and mRNAs, we calculated the Pearson correlation coefficients and the top 10 most relevant mRNAs were considered to be the co-expression partners of each lncRNA. The visible picture of co-expression network about paired lncRNAs and mRNAs was drawn using igraph program of R software. To predict the metabolic pathways and potential functions of lncRNAs, we performed the functional enrichment analysis for the co-expression mRNAs to find out significantly enriched Gene Ontology (GO) terms and Kyoto Ency- clopedia of Genes and Genomes (KEGG) pathway with the use of cluster Profiler software in R-version 4.0.4.

      Establish the Genomic Instability-Related Risk Score Model (GIRS)

      The clinical data of ccRCC patients and expression data of GILncRs were combined and all samples were divided into 2 sets randomly. One set was defined as training set to establish a prognostic model, the other one was defined as testing set to demonstrate the reliability of prognostic model. The bias of clinical parameters (such as gender, age, TNM stage) between 2 sets was examined by Chi-square test.
      GILncRs related to clinical outcomes of ccRCC patients were considered as survival-related GILncRs (sGILncRs). Univariate and multivariate Cox proportional hazard regression analysis were used to identify sGILncRs (P < .05) and measure the relationship between the expression level of sGILncRs and OS. Composed of the expression level of sGILncRs and the coefficients from the multivariate Cox regression analysis, a genomic instability-related risk score (GIRS) for prognosis prediction was constructed as follows:
      GIRS=i=1ncoef(lncRNAi)*exp(lncRNAi)


      In the formula, GIRS represents a risk score for ccRCC patients to predict the clinical outcome. LncRNAi is the ith sGILncR. Coef (lncRNAi) obtained from the multivariate Cox regression coefficient represents the contribution of lncRNAi to GIRS. Exp (lncRNAi) represents the expression level of lncRNAi . In the training set, the median GIRS of all patients was chosen as a critical point to divide patients into the low-risk group with low GIRS or high-risk group with high GIRS. To further explore the relationship between GIRS and clinicopathological characteristics of ccRCC, various clinical subgroups of ccRCC patients were analyzed such as gender, age, grade, neoadjuvant therapy and TNM stage.

      Kaplan-Meier Survival Estimation and ROC Curve

      The median survival and survival rate for each clinical subgroup were measured with the Kaplan-Meier analysis, and the difference of OS between the high-risk group and the low-risk group was evaluated by the log-rank test with a significant level of 5%. To verify the independence of GIRS from other clinical parameters, we used stratified analysis and multivariate Cox regression analysis to calculate Hazard ratio (HR) and 95% confidence interval (CI). We evaluated the performance of GIRS by the survival ROC curve using the survival ROC package in R-version 4.0.4, and the value of area under curve (AUC) was calculated. The prognostic correlation of sGILncRs was also confirmed in the GEO database. In the receiver operating characteristic (ROC) curve, the X-axis represents the false positive rate and the Y-axis represents the true positive rate.
      In addition, the GIRS model was compared with previous prognostic models established by Elena and Liu through ROC curve.

      Results

      Verification of Genomic Instability-Related lncRNAs (GILncRs) in Clear Cell Renal Cell Cancer Patients

      RNA-sequencing expression data, clinical files and mutation information were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/).
      • Ricketts C.J.
      • De Cubas A.A.
      • Fan H.
      • et al.
      The cancer genome atlas comprehensive molecular characterization of renal cell carcinoma.
      We distinguished mRNAs and lncRNAs through analyzing RNA-sequencing expression data. Subsequently, the mutation frequency of each sample was counted, and was sorted in ascending order. According to the mutation frequency, the top 25% of patients were classified as low mutation (LM) group, and the last 25% were classified as high mutation (HM) group. Hence, we can divide all the samples into the above 2 groups. We merged the lncRNA expression files and the sample mutation frequency files to obtain the expression level of lncRNAs in HM group and LM group, and then we try to find significantly differentially expressed lncRNAs which were defined as GILncRs between the 2 groups taking |logFC| > 1 and FDR-adjusted P-value < .05 as a standard (Supplementary Table 1). The heat map where the x axis represents the samples, and the y axis indicates the differential lncRNAs was drawn to show the difference of GILncRs expression between genomic unstable (GU) group and genomic stable (GS) group. As shown in Supplementary Figure 1, LINC00284, ZNF350-AS1, AL121820.2 and the other 21 lncRNAs are highly expressed in LM group while AC079466.1, AL035446.1, LINC00460 and the other 13 lncRNAs are highly expressed in HM group.
      We performed hierarchical clustering analysis using 539 samples from TCGA depending on the expression level of GILncRs. As we can see from Figure 1A, we categorized 539 samples into 2 groups based on the significantly different somatic mutation counts. The group which contained higher somatic mutation counts was defined as GU group while the other was defined as GS group. As shown in Figure 1B, there did exist notably significance in the median value of somatic mutation count between the GU group and the GS group (P = 4.2 e−11). As we can see, the GU group had higher somatic mutation count. Subsequently, we compared the expression level of PBRM1 gene (the second most highly mutated tumor suppressor gene in kidney cancer) between the GS group and GU group.
      • Cai W.
      • Su L.
      • Liao L.
      • et al.
      PBRM1 acts as a p53 lysine-acetylation reader to suppress renal tumor growth.
      ,
      • Espana-Agusti J.
      • Warren A.
      • Chew S.K.
      • Adams D.J.
      • Matakidou A.
      Loss of PBRM1 rescues VHL dependent replication stress to promote renal carcinogenesis.
      From Figure 1C, we can make a conclusion that the tumor suppressor gene PBRM1, which is specific for kidney cancer, is highly expressed in the GS group and is low expressed in the GU group. The 10 mRNAs that are most significantly related with GILncRs (P < .05) were used as target genes to draw a lncRNAs-mRNA co-expression network (Figure 2A). Supplementary Table 2 shows that AP000757.1, one example of GILncRs, is most correlated with MUC20, FAM189A2, ERP27, DHRS7, RHBG, NTN1, FOXI1, PIP, TBC1D1, ATP6V0D2 mRNAs. Gene Ontology (GO) analysis of GILncR-correlated mRNAs uncovered that these mRNAs play an important role in cellular monovalent inorganic cation homeostasis (Figure 2B), while Kyoto Ency- clopedia of Genes and Genomes (KEGG) pathway analysis revealed that mRNAs are related with collecting duct acid secretion and synaptic vesicle cycle (Figure 2C).
      Figure 1
      Figure 1Clustering of clear cell renal cell cancer (ccRCC) patients based on the expression level of GILncRs. (A) The blue and red cluster respectively represent genomic stable (GS) group and genomic unstable (GU) group. (B) Boxplots of somatic mutation counts in the two clusters, and the GU group has a significantly higher somatic mutation count than the GS group. (C) Boxplots of PBRM1 expression level in two clusters, and the GU group has a significantly lower expression level of PBRM1 than the GS group.
      Figure 2
      Figure 2Functional annotations of GILncRs. (A) Visual network of GILncRs and co-expressed mRNAs, where the blue circles are lncRNAs and the red circles are mRNAs. The co-expressed mRNA partners were performed by functional enrichment analysis of Kyoto Ency- clopedia of Genes and Genomes (KEGG) (B) and Gene Ontology (GO) (C)

      Establishment of Genomic Instability-Related Risk Score (GIRS) Model for Prediction of Clinical Outcomes

      LncRNA expression file was merged with the clinical data containing survival time, and then all the samples were randomly divided into 2 sets, such as the training set and the test set. It can be seen from Table 1 that the training set and the test set show no difference in clinical characteristics such as age, gender, grade, treatment and TNM stage, which can verify the randomness of our grouping. We conducted univariate cox regression analysis to analyze whether there was a relationship between the relative expression of 508 GILncRs and overall survival (OS) in the training set. The lncRNAs were defined as high-risk lncRNAs if the hazard ratio (HR) value is greater than 1, otherwise as low-risk lncRNAs. As shown in Supplementary Figure 2, the 6 lncRNAs are both high-risk lncRNAs. Subsequently, multivariate cox regression analysis was performed to find out the lncRNAs with the independent prognostic value in the training set. Then the 3 lncRNAs were found out and were used to establish a genomic instability-related risk score (GIRS) model. By calculating the risk score through the GIRS model, we can obtain the hazard ratio (HR) values of all samples both in the training set and the test set. The formula of GIRS model: HR value = (0.127 × relative expression of LINC00460) + (0.140 × relative expression of AC156455.1) + (0.037 × relative expression of LINC01606). As shown in Table 2, we can conclude that the lncRNAs may be risky factors for their high expressions were associated with poor clinical outcomes since the coefficients of lncRNAs are positive in the formula. Comparing all HR values with the median HR value in the training set, we can estimate the risk levels of all the samples. If the HR value was higher than the median HR value of the training set, the sample was categorized into the high-risk group, otherwise was categorized into the low-risk group.
      Table 1Clinical Information for Three Sets of ccRCC Patients in This Study
      All SetTesting SetTraining Set
      Covariates(n=507)(n=252)(n=255)P-value
      Age,no(%).51
       < = 65336(66.27)163(64.68)173(67.84)
       >65171(33.73)89(35.32)82(32.16)
      Gender, no(%).573
       FEMALE174(34.32)90(35.71)84(32.94)
       MALE333(65.68)162(64.29)171(67.06)
      Grade, no(%).391
       G1-2227(44.77)108(42.86)119(46.67)
       G3-4272(53.65)141(55.95)131(51.37)
       Unknow8(1.58)3(1.19)5(1.96)
      Stage, no(%).701
       Stage I-II306(60.36)155(61.51)151(59.22)
       Stage III-IV198(39.05)96(38.10)102(40.00)
       Unknow3(0.59)1(0.40)2(0.78)
      T,no(%).168
       T1-2324(63.91)169(67.06)155(60.78)
       T3-4183(36.09)83(32.94)100(39.22)
      M,no(%).818
       M0401(79.09)199(78.97)202(79.22)
       M178(15.38)37(14.68)41(16.08)
       Unknow28(5.52)16(6.35)12(4.71)
      N,no(%).254
       N0225(44.38)125(49.60)100(39.22)
       N1-316(3.16)6(2.38)10(3.92)
       Unknow266(52.47)121(48.02)145(56.86)
      Treatment, no(%).818
       Without neoadjuvant therapy491(96.84%)245(97.22%)246(96.47%)
       With neoadjuvant therapy16(3.16%)7(2.78%)9(3.53%)
      Table 2The Survival-Related GILncRs (sGILncRs) by Multivariate Cox Regression Analysis of GIRS Model
      IDCoefficientHRHR.95LHR.95HP-value
      LINC004600.1271.1361.0751.199<.001
      AC156455.10.1401.1511.0661.242<.001
      LINC016060.0371.0381.0071.070.017
      What's more, we performed Kaplan-Meier analysis to see if there are differences between the OS in the high-risk and low-risk groups. As shown in Figure 3A,B,C, the high-risk group was marked in red and the low-risk group was in blue, and we can make a conclusion that patients with lower risk have extremely better clinical outcomes. The conclusion reflects that the GIRS model can predict the prognosis of patients by means of dividing the patients into the above 2 groups. Additionally, the receiver operating characteristic (ROC) curve analysis of the GIRS model yielded an area under the curve (AUC). In our model, AUC was 0.686, 0.686, 0.687 in the training set, the test set and all samples respectively (Figure 3D,E,F), indicating that applying GIRS model to predict the patient’s survival is comparatively reliable.
      Figure 3
      Figure 3Effectiveness evaluations of genomic instability-related risk score (GIRS) model. Kaplan–Meier estimates of patients’ overall survival with low or high risk predicted by the GIRS model in training set (A), testing set (B) and all set (C). The receiver operating characteristic (ROC) curves of training set (D), testing set (E) and all set (F).

      Relationship between the Risk Score and Genomic Instability or lncRNA Expression

      With the purpose of figuring out the relationship between the risk and lncRNA expression, as well as genomic instability, several graphs were drawn to observe how the expression of lncRNAs in the GIRS model, the expression of PBRM1 and the somatic mutation count change with increasing risk scores. As shown in Figure 4A, the expression level of the high-risk lncRNA LINC00460 was upregulated with the increasing risk score, which indicated that LINC00460 is a high-risk lncRNA. Moreover, lncRNA AC156455.1 showed the same expression patterns while lncRNA LINC01606 remained unclear. As shown in Figure 4B, PBRM1 revealed lower expression level in the high-risk patients while the patients with lower risk expressed higher extent of PBRM1 since PBRM1 represented as a tumor suppressor gene (P = 2.2e−6). As we can see from Figure 4C, patients with higher risk were found to have notably higher somatic mutation count compared with those with lower risk (P = .035). In a word, patients with different extent of risk showed different patterns in symbols of genomic instability such as somatic mutation count and PBRM1 expression.
      Figure 4
      Figure 4The relationship between risk score model and genomic instability. (A) sGILncR expression patterns, somatic mutation count and expression of PBRM1 with increasing risk score. The expression of PBRM1 (B) and somatic mutation count (C) in the high-risk group and low-risk group of the ccRCC patients. The red and blue color respectively represent high-risk group and low-risk group. (D) The ROC curves of overall survival for the GIRS model, lncsig1 and lncsig2.

      Superiority of the GIRS Model to Published lncRNA Signatures in Survival Prediction

      Subsequently, we compared the GIRS model with other published lncRNA signatures, such as lncsig1 derived from Elena's study
      • Martens-Uzunova E.S.
      • Bottcher R.
      • Croce C.M.
      • Jenster G.
      • Visakorpi T.
      • Calin G.A.
      Long noncoding RNA in prostate, bladder, and kidney cancer.
      containing GAS5, H19, HIF1A-AS1, KCNQ1OT1, MALAT1, MEG3 lncRNAs and lncsig2 derived from Liu's study
      • Liu H.
      • Ye T.
      • Yang X.
      • et al.
      A panel of four-lncRNA signature as a potential biomarker for predicting survival in clear cell renal cell carcinoma.
      including TCL6, PVT1, MIR155HG, HAR1B lncRNAs. From the time-dependent ROC curves (Figure 4D), we can see that AUC of our model is larger than that of the other 2 signatures, illustrating that the GIRS model in survival prediction is more convincing.

      Independence of the GIRS Model from Other Clinical Features

      With the purpose of investigating whether our GIRS model was of prognostic value and was independent from other clinical factors, clinical features like age, gender, pathologic grade, pathologic stage along with our GIRS model were entered into multivariate cox regression analyses. We can see from Table 3, the GIRS model was highly correlated with OS when adjusted for age, gender, treatment, pathologic grade and pathologic stage (P < .001). Besides GIRS model, we found the other 3 clinical features (age, pathologic grade and pathologic stage) which were also observed to be of prognostic value in each set using the multivariate cox regression analysis. Hereafter, the patients were classified into various clinical groups, such as < = 65 year-old group (n = 336) and > 65 year-old group (n = 171), male group (n = 333) and female group (n = 174) (Supplementary Figure 3), T1-2 group (n = 324) and T3-4 group (n = 183), M0 group (n = 401) and M1 group (n = 78) (Supplementary Figure 4), stage Ⅰ/Ⅱ group (n = 306) and stage Ⅲ/Ⅳ group (n = 198), G1-2 group (n = 227) and G3-4 group (n = 272), With neoadjuvant therapy group (n = 16) and Without neoadjuvant therapy (n = 491) (Figure 5). Further analyses were conducted on the 2 separate risk groups to determine whether significant differences existed in OS in each clinical subset. As shown in Figure 5, our GIRS model was capable of serving as a prognostic factor independent from other clinical features in clear cell renal cell cancer, for the reason that there did notably differences exist in OS between the 2 separate risk group (P < .05).
      Table 3Univariate and Multivariate Cox Regression Analysis of GIRS Model and Overall Survival in Three ccRCC Patient Sets
      VariablesUnivariable ModelMultivariable Model
      HR95%CIP-valueHR95%CIP-value
      All set (n=507)age1.0291.015-1.043<.0011.0311.016-1.046<.001
      gender0.9690.702-1.338.849
      grade2.3081.864-2.858<.0011.4711.162-1.862.001
      stage1.9411.693-2.224<.0011.7181.472-2.006<.001
      treatment1.9630.969-3.849.053
      riskScore1.0261.016-1.037<.0011.0191.009-1.029<.001
      Training set (n=255)age1.0170.998-1.037.085
      gender0.8990.558-1.447.660
      grade2.3801.770-3.200<.0011.5631.110-2.201.010
      stage1.9041.561-2.323<.0011.6011.269-2.020<.001
      treatment2.0090.810-4.984.132
      riskScore1.0431.024-1.063<.0011.0311.010-1.051.003
      Testing set (n=252)age1.0411.021-1.061<.0011.0511.028-1.075<.001
      gender1.0440.672-1.624.847
      grade2.2131.623-3.019<.0011.4030.996-1.977.053
      stage1.9891.645-2.403<.0011.9041.539-2.355<.001
      treatment1.8700.683-5.115.223
      riskScore1.0211.010-1.031<.0011.0131.003-1.024.010
      Figure 5
      Figure 5Kaplan-Meier survival curve of subgroups of ccRCC patients with different grade (A, B), stage (C, D) and therapy (E, F). The high-risk group has worse outcome commonly except (F).

      External Prognostic Validation of the GIRS Model in GEO Database

      With the purpose of confirming the prognostic value of the GIRS model, we examined the association of lncRNAs in the GIRS model with clinical features in the independent dataset in The Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). We found that only one lncRNA (LINC00460) of 3 lncRNAs in the GIRS model was covered by GSE22541 which had large sample size and common clinicopathological characteristics. As shown in Figure 6A, the expression of LINC00460 was correlated with T staging and it is apparent that T3-4 stage revealed higher expression of LINC00460. However, differential expression of LINC00460 among other clinical subgroups were not reflected (Figure 6B, C and D). We tried to find the other 2 lncRNAs, AC156455.1 and LINC01606 in the datasets of GEO database, but failed. Moreover, we divided the samples in the dataset GSE22541 into 2 groups based on the average expression level of LINC00460 and Kaplan-Meier analysis was performed to see whether there existed difference in OS between the high-expression group and the low-expression group. It can be seen that the patients with higher expression of LINC00460 had a poorer prognosis than those with lower expression although the difference was not significant enough (P = .104, Figure 6E). We suppose that the biased results attributed to the small sample size, hence we plotted the survival curve of LINC00460 on the Kaplan-Meier Plotter (https://kmplot.com/analysis/). As shown in Figure 6F, patients with different LINC00460 expression tend to have distinct clinical outcomes showing that LINC00460 exhibited significant prognostic potential.
      Figure 6
      Figure 6Association of lncRNAs in the GIRS model with clinical features characteristics in The Gene Expression Omnibus (GEO) database. The expression of LINC00460 in different T-stage (A), N-stage (B), M-stage (C) and grade (D) within GEO data. (E)Survival curve of ccRCC patients in GEO database. Kaplan-Meier analysis of two group based on the expression level of LINC00460 by R-version 4.0.4. (F) Survival curve drawn on the Kaplan-Meier Plotter. The red line and blue line respectively represent the high-expression group and low-expression group.

      Discussion

      Although improvements in diagnosis and treatment of renal cell cancer (RCC) have occurred in the last 20 years, renal cell cancer continues to have significant morbidity and mortality. About a quarter of RCC patients are diagnosed as late-stage patients where metastasis of tumor can often been observed. The 5-year survival rate of patients with stage I, II, III and IV renal cancer after treatment reached 92%, 86%, 64% and 23%, respectively. The general principle of treatment for RCC suggests that surgical treatment ought to be adopted for early-stage patients, while systematic treatment should be based on internal medicine for the patients with advanced cancer. Since renal cell cancer is naturally insensitive to radiotherapy and chemotherapy, biological targeted therapy or immunotherapy is generally used.
      • Kalimutho M.
      • Nones K.
      • Srihari S.
      • Duijf P.H.G.
      • Waddell N.
      • Khanna K.K.
      Patterns of genomic instability in breast cancer.
      There are no specific predictors for us to decide whether targeted therapy should be adopted. Presently, we adopt biological targeted therapy according to pathological type, MSKCC score, drug adverse reaction and patient's will. Besides, markers for predicting the efficacy of targeted therapy still remain unclear.
      Genomic instability (GI) represents a notably indispensable feature of cancer. However, genomic instability is hard to be measured quantitatively, so people are committed to identifying small RNAs associated with genomic instability and developing gene or small RNA models to predict genomic instability. Recently, a class of novel non-coding ribonucleotides (lncRNAs) were known to be important components of tumor biology. Emerging evidence reveals that lncRNAs play a crucial role in maintain genome stability and specific lncRNAs are essential for survival of tumor cell and prevention of tumorigenesis, which reminds us that lncRNAs may become a valuable tool for clinical prognosis and therapeutic decision-making.
      Then, we combined the lncRNA expression profile of ccRCC with somatic mutation file to obtain the expression of lncRNA in high mutation (HM) group and low mutation (LM) group, and analyzed the difference. The 47 lncRNAs with significant difference between the above 2 groups were screened out, indicating that the 47 lncRNAs were correlated with genomic instability. We further investigated whether GILncRs can predict clinical outcomes, therefore we established the genomic instability-related risk score (GIRS) model containing 5 genomic instability-related lncRNAs. In the training set, GIRS model divided the patients into 2 separate risk groups and we analyzed that there were significant differences in survival rates between the above 2 groups, which was also observed in the independent test set. Moreover, GIRS was significantly associated with PBRM1 expression in clear cell renal cell cancer, while PBRM1 is an important indicator of genomic instability. Through literature research, we found the biological functions of 2 lncRNAs (LINC01606, LINC00460) in the GIRS model. LINC01606 which is located at chromosome8 (chromosome8:57218276-57232954) tend to have a relationship with metastasis and invasion of gastric cancer as previously reported.
      • Luo Y.
      • Tan W.
      • Jia W.
      • et al.
      The long non-coding RNA LINC01606 contributes to the metastasis and invasion of human gastric cancer and is associated with Wnt/beta-catenin signaling.
      Additionally, LINC00460 has a strong correlation with poor prognosis of patients. Moreover, our GIRS model were entered into multivariate cox regression analysis with a constant range of P value (less than .05), indicating that our GIRS model was capable of serving as a prognostic factor independent from other clinical features in ccRCC. What's more, it was apparent that significant differences existed in OS in each clinical subset, revealing that our model is suitable for patients of any sex, any treatment, all ages, all stages, all pathological grades and almost all TNM stages. However, we didn't observe differences in OS in patients with neoadjuvant therapy due to the insufficient sample size. Therefore, the above findings are essential for ongoing research in precision medicine and immunonomonology, helping identifying patients who are more likely to suffer from ccRCC but sensitive to targeted therapy or immunotherapy. It also provides basic information for the development of advisable combination therapy to overcome drug resistance. Although our study endeavors to use genomic instability and lncRNAs to predict prognosis, there still exist some limitations that require further study.eg, it is confusing that the difference in somatic mutation count between the 2 separate risk groups does not show extremely significantly but extraordinarily significant difference exists in the PBRM1 expression of the 2 risk groups, which can be explained by the comment of Jennifer Couzin-Frankel. The researcher believes that some mutations may have more effects than others although cancer occurs when various mutations accumulate to a certain extent.
      • Couzin-Frankel J.
      Science communication. backlash greets 'bad luck' cancer study and coverage.
      In our study, it is apparent that PBRM1 mutation instead of the somatic mutation constitutes a more essential part in the occurrence of clear cell renal cell cancer. What's more, we only found LINC00460 in the GEO dataset and then we analyzed the correlation of LINC00460 and overall survival, while the biological function of lncRNAs in the GIRS model still remained to be further explored. Furthermore, we need more independent data sets to validate GIRS to ensure its robustness and repeatability.

      Conclusions

      In conclusion, we verified the effects of genomic instability related lncRNAs on predicting the clinical outcomes of patients. We rest hope on the results which may be capable of developing a credible and referable risk assessing model and making remarkable effort to provide new insight into the targeted therapy of clear cell renal cancer.
      • Sanchez-Gastaldo A.
      • Kempf E.
      • Gonzalez Del Alba A.
      • Duran I.
      Systemic treatment of renal cell cancer: A comprehensive review.
      ,
      • Chen W.
      • Hill H.
      • Christie A.
      • et al.
      Targeting renal cell carcinoma with a HIF-2 antagonist.
      • Courtney K.D.
      • Ma Y.
      • Diaz de Leon A.
      • et al.
      HIF-2 complex dissociation, target inhibition, and acquired resistance with PT2385, a first-in-class HIF-2 inhibitor, in patients with clear cell renal cell Carcinoma.
      • Miao D.
      • Margolis C.A.
      • Gao W.
      • et al.
      Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma.

      Author Contributions

      S.L. and H.H. conceived and designed the research; S.L., H.H. and Z.S. collected and analyzed the data; H.H. and C.D. made figures and tables; H.H. and S.L. contributed to the writing of the first draft; L.Q. and L.Y. supervised the implementation of the project and reviewed the manuscript. All authors have read and approved the submitted version of the manuscript.

      Institutional Review Board Statement

      This study is based on open-source data, so there are no ethical issues to be reviewed.

      Informed Consent Statement

      This study is based on TCGA and GEO that belong to public databases. The patients involved in the database have obtained ethical approval already.

      Data Availability Statement

      The RNA-sequencing data and clinical information of ccRCC patients are available in TCGA and GEO with accession number GSE22541.

      Acknowledgments

      We express our best gratitude to TCGA and GEO databases for providing their platforms and contributors for uploading their precious datasets. We are grateful to Yuanyuan Zhang for the helpful technical guidance on data processing.

      Disclosure

      The authors declare no conflict of interest in this study.

      Appendix. Supplementary materials

      References

        • Capitanio U.
        • Bensalah K.
        • Bex A.
        • et al.
        Epidemiology of renal cell carcinoma.
        Eur Urol. 2019; 75: 74-84https://doi.org/10.1016/j.eururo.2018.08.036
        • O'Rourke C.J.
        • Knabben V.
        • Bolton E.
        • et al.
        Manipulating the epigenome for the treatment of urological malignancies.
        Pharmacol Ther. 2013; 138: 185-196https://doi.org/10.1016/j.pharmthera.2013.01.007
        • Linehan W.M.
        Genetic basis of kidney cancer: role of genomics for the development of disease-based therapeutics.
        Genome Res. 2012; 22: 2089-2100https://doi.org/10.1101/gr.131110.111
        • Vuong L.
        • Kotecha R.R.
        • Voss M.H.
        • Hakimi A.A.
        Tumor microenvironment dynamics in clear-cell renal cell carcinoma.
        Cancer Discov. 2019; 9: 1349-1357https://doi.org/10.1158/2159-8290.CD-19-0499
        • Rydzanicz M.
        • Wrzesinski T.
        • Bluyssen H.A.
        • Wesoly J.
        Genomics and epigenomics of clear cell renal cell carcinoma: recent developments and potential applications.
        Cancer Lett. 2013; 341: 111-126https://doi.org/10.1016/j.canlet.2013.08.006
        • Song M.
        Recent developments in small molecule therapies for renal cell carcinoma.
        Eur J Med Chem. 2017; 142: 383-392https://doi.org/10.1016/j.ejmech.2017.08.007
        • Sanchez-Gastaldo A.
        • Kempf E.
        • Gonzalez Del Alba A.
        • Duran I.
        Systemic treatment of renal cell cancer: A comprehensive review.
        Cancer Treat Rev. 2017; 60: 77-89https://doi.org/10.1016/j.ctrv.2017.08.010
        • Guo C.
        • Zhao H.
        • Wang Y.
        • et al.
        Prognostic value of the neo-immunoscore in renal cell carcinoma.
        Front Oncol. 2019; 9: 439https://doi.org/10.3389/fonc.2019.00439
        • Negrini S.
        • Gorgoulis V.G.
        • Halazonetis T.D.
        Genomic instability–an evolving hallmark of cancer.
        Nat Rev Mol Cell Biol. 2010; 11: 220-228https://doi.org/10.1038/nrm2858
        • Kalimutho M.
        • Nones K.
        • Srihari S.
        • Duijf P.H.G.
        • Waddell N.
        • Khanna K.K.
        Patterns of genomic instability in breast cancer.
        Trends Pharmacol Sci. 2019; 40: 198-211https://doi.org/10.1016/j.tips.2019.01.005
        • Poole L.A.
        • Cortez D.
        Functions of SMARCAL1, ZRANB3, and HLTF in maintaining genome stability.
        Crit Rev Biochem Mol Biol. 2017; 52: 696-714https://doi.org/10.1080/10409238.2017.1380597
        • Syed A.
        • Tainer J.A.
        The MRE11-RAD50-NBS1 complex conducts the orchestration of damage signaling and outcomes to stress in DNA replication and repair.
        Annu Rev Biochem. 2018; 87: 263-294https://doi.org/10.1146/annurev-biochem-062917-012415
        • Cheetham S.W.
        • Gruhl F.
        • Mattick J.S.
        • Dinger M.E.
        Long noncoding RNAs and the genetics of cancer.
        Br J Cancer. 2013; 108: 2419-2425https://doi.org/10.1038/bjc.2013.233
        • Moghaddas Sani H.
        • Hejazian M.
        • Hosseinian Khatibi S.M.
        • Ardalan M.
        • Zununi Vahed S.
        Long non-coding RNAs: An essential emerging field in kidney pathogenesis.
        Biomed Pharmacother. 2018; 99: 755-765https://doi.org/10.1016/j.biopha.2018.01.122
        • Martens-Uzunova E.S.
        • Bottcher R.
        • Croce C.M.
        • Jenster G.
        • Visakorpi T.
        • Calin G.A.
        Long noncoding RNA in prostate, bladder, and kidney cancer.
        Eur Urol. 2014; 65: 1140-1151https://doi.org/10.1016/j.eururo.2013.12.003
        • Lee S.
        • Kopp F.
        • Chang T.C.
        • et al.
        Noncoding RNA NORAD regulates genomic stability by sequestering PUMILIO proteins.
        Cell. 2016; 164: 69-80https://doi.org/10.1016/j.cell.2015.12.017
        • Munschauer M.
        • Nguyen C.T.
        • Sirokman K.
        • et al.
        The NORAD lncRNA assembles a topoisomerase complex critical for genome stability.
        Nature. 2018; 561: 132-136https://doi.org/10.1038/s41586-018-0453-z
        • Hu W.L.
        • Jin L.
        • Xu A.
        • et al.
        GUARDIN is a p53-responsive long non-coding RNA that is essential for genomic stability.
        Nat Cell Biol. 2018; 20: 492-502https://doi.org/10.1038/s41556-018-0066-7
        • Khanduja J.S.
        • Calvo I.A.
        • Joh R.I.
        • Hill I.T.
        • Motamedi M.
        Nuclear noncoding RNAs and genome stability.
        Mol Cell. 2016; 63: 7-20https://doi.org/10.1016/j.molcel.2016.06.011
        • Khorkova O.
        • Hsiao J.
        • Wahlestedt C.
        Basic biology and therapeutic implications of lncRNA.
        Adv Drug Deliv Rev. 2015; 87: 15-24https://doi.org/10.1016/j.addr.2015.05.012
        • Rothschild G.
        • Basu U.
        Lingering questions about enhancer RNA and enhancer transcription-coupled genomic instability.
        Trends Genet. 2017; 33: 143-154https://doi.org/10.1016/j.tig.2016.12.002
        • Wan G.
        • Liu Y.
        • Han C.
        • Zhang X.
        • Lu X.
        Noncoding RNAs in DNA repair and genome integrity.
        Antioxid Redox Signal. 2014; 20: 655-677https://doi.org/10.1089/ars.2013.5514
        • Ricketts C.J.
        • De Cubas A.A.
        • Fan H.
        • et al.
        The cancer genome atlas comprehensive molecular characterization of renal cell carcinoma.
        Cell Rep. 2018; 23 (e315): 313-326https://doi.org/10.1016/j.celrep.2018.03.075
        • Cai W.
        • Su L.
        • Liao L.
        • et al.
        PBRM1 acts as a p53 lysine-acetylation reader to suppress renal tumor growth.
        Nat Commun. 2019; 10: 5800https://doi.org/10.1038/s41467-019-13608-1
        • Espana-Agusti J.
        • Warren A.
        • Chew S.K.
        • Adams D.J.
        • Matakidou A.
        Loss of PBRM1 rescues VHL dependent replication stress to promote renal carcinogenesis.
        Nat Commun. 2017; 8: 2026https://doi.org/10.1038/s41467-017-02245-1
        • Liu H.
        • Ye T.
        • Yang X.
        • et al.
        A panel of four-lncRNA signature as a potential biomarker for predicting survival in clear cell renal cell carcinoma.
        J Cancer. 2020; 11: 4274-4283https://doi.org/10.7150/jca.40421
        • Luo Y.
        • Tan W.
        • Jia W.
        • et al.
        The long non-coding RNA LINC01606 contributes to the metastasis and invasion of human gastric cancer and is associated with Wnt/beta-catenin signaling.
        Int J Biochem Cell Biol. 2018; 103: 125-134https://doi.org/10.1016/j.biocel.2018.08.012
        • Couzin-Frankel J.
        Science communication. backlash greets 'bad luck' cancer study and coverage.
        Science. 2015; 347: 224https://doi.org/10.1126/science.347.6219.224
        • Chen W.
        • Hill H.
        • Christie A.
        • et al.
        Targeting renal cell carcinoma with a HIF-2 antagonist.
        Nature. 2016; 539: 112-117https://doi.org/10.1038/nature19796
        • Courtney K.D.
        • Ma Y.
        • Diaz de Leon A.
        • et al.
        HIF-2 complex dissociation, target inhibition, and acquired resistance with PT2385, a first-in-class HIF-2 inhibitor, in patients with clear cell renal cell Carcinoma.
        Clin Cancer Res. 2020; 26: 793-803https://doi.org/10.1158/1078-0432.CCR-19-1459
        • Miao D.
        • Margolis C.A.
        • Gao W.
        • et al.
        Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma.
        Science. 2018; 359: 801-806https://doi.org/10.1126/science.aan5951