The dataset for KIRC was downloaded from The Cancer Genome Atlas (TCGA) located at UCSC Xena (https://xenabrowser.net/). RNA sequencing (RNA-seq) data was extracted from the TCGA data portal. Values for fragments per kilobase million (FPKM) were converted to transcripts per million (TPM).

A total of 130 purine-related genes were obtained from KEGG (hsa00230; purine metabolism), and 129 overlapping genes were identified by intersecting with the TCGA gene set. Single-factor Cox filtering was used to identify 65 genes with p 0.01. NMF clustering was used to classify KIRC and to identify purine-related patterns, with the patients then grouped for subsequent analysis.

The limma package (version 3.56.2, https://bioinf.wehi.edu.au/limma/) was employed for single-factor filtering and to identify differentially expressed lncRNAs in the purine-associated category, using a significance threshold of p 0.05 and $|$logFC$|$ $>$1. Random survival forest was then used to perform further screening, and the Lasso method was employed to build a model using the selected genes. The risk score was calculated by multiplying gene expression values with the regression coefficients. The surv_cutpoint function of the survminer package (version 0.4.9, https://rpkgs.datanovia.com/survminer/index.html) was then applied to classify patients into high- and low-risk groups. The timeROC R package (version 0.4, https://www.rdocumentation.org/packages/timeROC/versions/0.4) was used to plot time-ROC curves.

The MCPcounter and TIMER algorithms were utilized to evaluate immune cell abundance in KIRC samples and to identify disparities in immune cell infiltration across distinct clustering categories or risk groups. Gene Set Enrichment Analysis (GSEA) was performed to examine KEGG pathway regulation. The maftools R package (version 2.16.0, https://github.com/PoisonAlien/maftools) was used for mutation analysis.

The cells of the human RCC line 786-O (AW-CCH060) and CAKI-1 (AW-CCH173) were purchased from Abiowell (Changsha, Hunan, China). 786-O or CAKI-1 cells were cultured in RPMI-1640 or McCoy’s 5A medium containing 10% fetal bovine serum and 1% Penicillin/Streptomycin. All cells were maintained at 37 °C with 5% CO${}_2$ in a humidified atmosphere. All cell lines have been authenticated in the past three years. All cell lines were identified by STR profile. All experiments were performed with mycoplasma-free cells. Both knockdown and overexpression of LINC01671 were performed. 786-O cells were divided into control, sh-NC, and sh-LINC01671 groups, and CAKI-1 cells into control, oe-NC, and oe-LINC01671 groups. The sh-LINC01671, oe-LINC01671, and negative controls (sh-NC and oe-NC) were provided by HonroGene (Changsha, Hunan, China). All transfections were performed using Lipofectamine 2000 (11668019, Invitrogen, Waltham, MA, USA).

qRT-PCR was utilized to evaluate LINC01671, ARAP1-AS1, and LINC02747 levels. Total RNA was extracted and reverse transcribed into cDNAs. Ultra SYBR Mixture (CW2601, CWBIO, Cambridge, MA, USA) was used to test on the ABI 7900 system. Gene expression was calculated using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ method, with GAPDH as the internal reference. Primer sequences were: LINC01671-F: TCAGGAACACCTCACAGGTC, LINC01671-R: GCAAACTCCAAGAGGAGTCCA; ARAP1-AS1-F: TCCTCTACAGCACCCGCTTT, ARAP1-AS1-R: CCACCCTTTCAGAGGCGTGAG; LINC02747-F: GAAGATGTGCACCTGCCGAG, LINC02747-R: GGTTGAGTTCAATGGCAGCA; GAPDH-F: ACAGCCTCAAGATCATCAGC, GAPDH-R: GGTCATGAGTCCTTCCACGAT.

Cells were seeded at a density of 1 $\times$ 10${}^4$/100 µL in a 96-well plate and incubated at 37 ℃ in 5% CO${}_2$. After adding 10 µL CCK-8 (NU679, DOJINDO, Tokyo, Japan), cells were incubated at 37 ℃ in 5% CO${}_2$ for 4 h. Optical density (OD) values for absorbance at 450 nm were analyzed using an enzyme marker (MB-530, Heales, Shenzhen, Guangdong, China).

Cells were suspended in serum-free medium at 1 $\times$ 10${}^6$/mL and 100 µL was added to the upper chamber of a Transwell (33318035, Corning, Somerville, MA, USA). The lower chamber was filled with complete medium containing 10% fetal bovine serum (FBS). After removal of the culture medium from the upper chamber, the upper surface of cells was wiped with a wet cotton swab. The cells were then fixed and stained with crystal violet. Cells on the outer surface of the upper chamber were observed and photographed under a microscope (Olympus, Tokyo, Japan).

Cells were digested and centrifuged, and about 3.2 $\times$ 10${}^5$ cells were collected. These were suspended in 500 µL of binding buffer, and 5 µL allophycocyanin (APC) and propidium iodide were added and mixed. Reaction was carried out in the dark, and the cells analyzed by flow cytometry (A00-1-1102, Beckman, Pasadena, CA, USA) within 1 h.

A TUNEL apoptosis detection kit (FITC) (40306ES50, Yeasen, Shanghai, China) was used to evaluate cell apoptosis. After fixation, cells were permeabilized with Triton X-100 and sodium citrate solution, and then treated with fluorescent-labeled nucleotides (dUTP) and TdT. TUNEL-positive cells with green fluorescence were observed and quantified with a fluorescence microscope.

Normative variables were tested using the Shapiro-Wilk test, while normally distributed variables were compared with unpaired Student’s t-test. Non-normally distributed variables were compared using the Wilcoxon test. One-way analysis of variance (ANOVA) was used as a parametric method to compare multiple groups, and the Kruskal-Wallis test as a non-parametric method. For each dataset, patients were categorized by binary risk score, with the R package ggplot2 (version 3.4.3, https://ggplot2.tidyverse.org/) used to visualize data. The Benjamini-Hochberg method was used to analyze differential gene expression. Significant genes were identified via the conversion of p-values to false discovery rate (FDR). The Kaplan-Meier method was used to compare survival of different patient groups, with the logarithmic rank test used to assess whether differences were statistically significant. All survival curves were generated using the R package survminer, and all heatmaps using pheatmap. Statistical analysis was performed using R (version 3.6.1, https://www.r-project.org/), with statistical significance set at p 0.05 for two-sided tests.

A total of 130 purine-related genes (KEGG: hsa00230; purine metabolism) were identified from the literature. The intersection of these genes with TCGA resulted in 129 genes. Subsequently, 65 genes were selected by univariate Cox analysis (p 0.01, Table 1). KIRC were grouped into Cluster-1 and Cluster-2 according to purine genes. As the number of clusters increased, the cophenetic coefficient decreased (Fig. 1A). Survival analysis showed that Cluster-2 patients had better survival than Cluster-1 patients (Fig. 1B). A volcano plot was used to visualize differences in lncRNAs between the two clusters. The limma package identified 22 differentially expressed lncRNAs between Cluster-1 and Cluster-2 (Fig. 1C). A heatmap further visualized the expression of purine-related lncRNAs between the two clusters (Fig. 1D).

Fig. 1.

Clustering of kidney renal clear cell carcinoma (KIRC) based on purine genes. (A) Non-negative matrix factorization (NMF) clustering analysis. (B) Survival analysis for Cluster-1 and Cluster-2. (C) Differentially expressed long non-coding RNAs (lncRNAs) between Cluster-1 and Cluster-2 were visualized by a volcano plot. (D) Clustering heatmap showing the expression of purine-related lncRNAs in Cluster-1 and Cluster-2. M, Metastasis; N, Node; T, Tumor. ****p 0.0001.

Purine-related lncRNAs were first screened by univariate Cox analysis (Fig. 2A). This identified ARAP1-AS1 with an increased hazard ratio, and 15 lncRNAs with a decreased hazard ratio (LINC01320, LINC02274, LINC00671, LINC02532, ADAMTS9-AS1, LINC01697, C6orf223, LHFPL3-AS2, DRAIC, LINC01508, PRKARIB-AS2, LINC2747, LINC02754, LINC01671, and LINC01550). Next, random forest analysis was used to reduce the dimension, resulting in 5 lncRNAs (LINC01671, ARAP1-AS1, LINC02747, ADAMTS9-AS1, and LINC01697; Fig. 2B). Lasso analysis then identified 3 lncRNAs (LINC01671, ARAP1-AS1 and LINC02747; Fig. 2C). Finally, the Lasso method was used to obtain a risk score model comprised of 3 lncRNAs: –0.1406 $\times$LINC01671 + 0.0739 $\times$ARAP1-AS1 – 0.1592 $\times$LINC02747. A clustering heatmap was used to visualize expression of the 3 lncRNAs (LINC01671, ARAP1-AS1, and LINC02747; Fig. 2D).

Fig. 2.

Development of a risk score (Purine Score) based on differentially expressed, purine-related lncRNAs. (A) Univariate Cox analysis identified 16 purine-related lncRNAs. (B) Random forest analysis identified 5 lncRNAs. (C) Lasso analysis identified 3 lncRNAs. (D) Clustering heatmap used to visualize expression of 3 lncRNAs. ****p 0.0001.

According to the risk score model established following TCGA survival analysis, patients with high risk scores had worse prognosis (p 0.05, Fig. 3A). Moreover, based on receiver operating characteristic (ROC) analysis, the area under the curve (AUC) for Purine Score in TCGA were 0.688, 0.681, and 0.695 for 1-, 3-, and 5-year true-positive rates, respectively (Fig. 3B). We also investigated prognostic factors using univariate and multivariate Cox analyses. Age and metastasis (M) were found to be independent prognostic factors for KIRC patients (Fig. 3C). Finally, the Purine Score was evaluated according to various clinical features (Fig. 3D). Significant differences in the Purine Score were found according to gender, grade, stage, tumor (T), node (N) and status (p 0.05).

Fig. 3.

Prediction of outcome in KIRC patients using Purine Score. (A) Survival analysis based on the risk score. (B) Receiver operating characteristic (ROC) analysis. (C) Identification of prognostic factors using univariate and multivariate Cox analyses. (D) Analysis of Purine Score according to clinical features. TCGA, The Cancer Genome Atlas.

Next, we performed a correlation analysis between prognosis and immune cell infiltration. As shown in Fig. 4A, the Purine Score correlated with cellular immune response and cell components, as determined by the MCPcounter and TIMER algorithms. The survival model included abnormalities in cytotoxic lymphocytes, NK cells, myeloid dendritic cells, monocytic lineage, endothelial cells, neutrophils, B cells, T cells CD4, neutrophils, macrophages, and DCs. Patients with a low Purine Score showed significantly higher scores for ESTIMATE (p = 0.0008), immune (p = 0.00026), and stromal (p = 0.0016) compared to patients with a high Purine Score (Fig. 4B). Correlation of immune regulation factors with the Purine Score are shown in Supplementary Fig. 1. The classification categories for immune checkpoints were cell adhesion, antigen presentation, co-stimulator, co-inhibitor, ligand, other, and receptor. A significant association was observed between the expression of immune checkpoint genes and the Purine Score.

Fig. 4.

Analysis of immune cell infiltration and immune checkpoints. (A) Analysis of immune cell infiltration. (B) ESTIMATE, Immune, and Stromal Scores. *p 0.05, ****p 0.0001.

Somatic mutation analysis revealed alterations in 77 of 97 (79.38%) KIRC with a high Purine Score (Missense Mutation, Frame Shift Del, Splice Site, Frame Shift Ins, In Frame Ins, Nonsense Mutation, Translation Start Site, and Multi Hit). Moreover, 211 of 235 (89.79%) KIRC with a low Purine Score showed alterations (Nonsense Mutation, Frame Shift Del, Frame Shift Ins, Missense Mutation, In Frame Ins, Translation Start Site, In Frame Del, and Nonstop Mutation, Fig. 5A). Gene mutation frequencies for the high and low Purine Score groups are shown in Fig. 5B. The mutation frequencies for TP53, TRIOBP, PKHD1, NPHP3, TLN2, CABIN1, ABCC6, XIRP2, and CHD4 were significantly higher in the high Purine Score group compared to the low Purine Score group, whereas the mutation frequencies for PBRM1 and VHL were significantly lower. Furthermore, in the high Purine Score group, VHL mutation co-occurred with PBRM1 mutation, PBRM1 with SETD2, SETD2 with MUC17, MTOR with CHD4, KDM5C with ABCC6, TP53 with MUC17, XIRP2 with CHD4 and CSMD3, CSMD3 with CABIN1, and CABIN1 with DST. In the low Purine Score group, PBRM1 mutation did not co-occur with BAP1, TTN mutation co-occurred with BRCA2, ANK3, HMCN1 and BAP1 mutation, SETD2 mutation co-occurred with LRP2, BAP1 and MUC16 mutation, and ARID1A with DNAH9 (Fig. 5C).

Fig. 5.

Analysis of gene mutations. (A) Somatic mutation analysis. (B) Gene mutation frequencies in the high and low Purine Score groups. (C) Co-occurrence of mutated genes in the high and low Purine Score groups. ◼p 0.05, *p 0.01.

Next, we used in vitro experiments to validate cell expression of the selected 3 lncRNAs (LINC01671, ARAP1-AS1, and LINC02274). qRT-PCR revealed high expression of LINC01671 in 786-O cells, and low expression in CAKI-1 cells. Since ARAP1-AS1 and LINC02274 were highly expressed in both 786-O and CAKI-1 cell lines (Fig. 6A), LINC01671 was therefore selected for further study. Survival analysis showed that high expression of LINC01671 was associated with improved survival (Fig. 6B). Next, LINC01671 was knocked down in 786-O cells and overexpressed in CAKI-1 cells. Successful transfection of sh-LINC01671 and of oe-LINC01671 was achieved, as shown in Fig. 6C. Cell function experiments showed that knockdown of LINC01671 in 786-O cells promoted their proliferation and migration, but inhibited apoptosis. In contrast, overexpression of LINC01671 in CAKI-1 cells inhibited their proliferation and migration, while promoting apoptosis (Fig. 6D–G).

Fig. 6.

Validation of lncRNA expression. (A) Quantitative real-time PCR (qRT-PCR) analysis of LINC01671, ARAP1-AS1, and LINC02274 levels in 786-O and CAKI-1 cells. (B) Survival analysis according to LINC01671 expression level. (C) qRT-PCR analysis of LINC01671 expression after transfection of 786-O and CAKI-1 cells. (D) Cell counting kit 8 (CCK-8) assay results for cell proliferation. (E) Transwell cell migration results. (F) Flow cytometry analysis of cell apoptosis. (G) Terminal Deoxynucleotidyl Transferase mediated dUTP Nick-End Labeling (TUNEL) analysis of cell apoptosis. * p 0.05 vs. sh-NC, # p 0.05 vs. oe-NC.

Finally, we performed functional enrichment analysis of LINC01671. GSEA showed that LINC01671 was mainly enriched in the MAPK (normalized enrichment score (NES) = 1.6, p 0.001), NF-kappa B (NES = 1.64, p 0.001), mTOR (NES = 1.56, p 0.001), PI3K-Akt (NES = 1.61, p 0.001) and Wnt (NES = 1.51, p = 0.0001) signaling pathways (Fig. 7). These results suggest that LINC01671 may positively regulate the MAPK, NF-kappa B, mTOR, PI3K-Akt and Wnt signaling pathways. The flow chart is shown in Supplementary Fig. 2.

Fig. 7.

Functional enrichment analysis of LINC01671. Gene Set Enrichment Analysis (GSEA) of the function and pathways for LINC01671.