Disulfideptosis-associated lncRNAs reveal features of prognostic, immune escape, tumor mutation, and tumor malignant progression in renal clear cell carcinoma

Acquisition of raw data

Clinical data of 533 renal clear cell carcinoma cases and expression profile data of all RNAs were downloaded from TCGA (The Cancer Genome Atlas). Similarly, 403 simple nucleotide variation data were downloaded from TCGA. Genes associated with Disulfideptosis from previous studies on Disulfideptosis [3], Tumor Immune Dysfunction and Exclusion (TIDE) data generated from the RNA transcriptome expression profiles of TCGA were downloaded from TIDE (http://tide.dfci.harvard.edu/) to calculate their TIDE values.

Obtaining Disulfideptosis-related lncRNAs

Based on the downloaded TCGA RNA transcriptome expression profile, we used Perl software (bictype.perl) to separate the RNA transcriptome expression profile into lncRNA and mRNA expression profiles. Then, using R software (packages “ggplot” and “ggalluvial”), we depicted the co-expression relationship between lncRNA and Disulfideptosis-associated genes via Pearson correlation analysis (correlation with p < 0.05) as a Sankey diagram.

Construction of the prognostic signature

Using R software (packages “survival”, “caret”, “glmnet”, “survminer”, and “timeROC”), we conducted univariate Cox regression analysis on Disulfideptosis-associated lncRNAs (significance criterion was p < 0.001), obtained Disulfideptosis-associated lncRNAs significantly associated with prognosis time, and plotted a forest map. Then we conducted a Lasso regression analysis (randomly dividing it into two groups, performing cross-validation, and finding the minimum error value). After that, we performed a multivariate Cox regression analysis of related factors, filtering out independent prognostic factors (5 disulfidptosis-related lncRNAs), thus constructing the prognostic signature. Using R software (packages “tidyverse”, “ggplot2”, “ggExtra”), we plotted a correlation heatmap of these 5 independent prognostic lncRNAs and 13 Disulfideptosis-associated genes through Pearson correlation analysis (correlation with p < 0.05).

Evaluation and validation of the risk model

Using R software (packages “survival”, “caret”, and “survminer”), we randomly divided all 533 samples into two groups (training group and testing group). We then conducted a chi-square test on the clinical characteristics (Age, Gender, Grade, and Stage) of both groups, comparing whether there is a significant difference in the clinical characteristics of the two groups. A p-value less than 0.05 is considered statistically significant. Using R software (package “pheatmap”), we plotted risk curves, survival status maps, and risk heatmaps for both the training and testing groups (divided into high and low-risk groups based on the median). Using R software (packages “survival” and “survminer”), we plotted Kaplan-Meier (KM) curves for high and low-risk groups in both the training and testing sets, including Overall Survival (OS) and Progression Free Survival (PFS). Using R software (the “survival” package), we conducted univariate and multivariate independent prognostic analyses for the risk model. Using R software (packages “survival”, “survminer”, and “timeROC”), we generated the 1, 3, and 5-year ROC curves for the risk model and the 5-year ROC curves for Age, Gender, Grade, and Stage, and made a comparison. Using R software (packages “limma” and “scatterplot3d”), we performed a Principal Component Analysis (PCA) on the risk model to see if it can distinguish all lncRNAs, Disulfideptosis-associated genes lncRNAs, Disulfideptosis-associated genes, and all genes and visualized the results.

Construction of the predictive nomogram

Using R software (packages “survival”, “regplot”, and “rms”), we incorporated the risk model into the construction of a predictive nomogram based on the risk score. Subsequently, we generated calibration curves for the 1, 3, and 5-year periods of the nomogram.

Clinical characteristics of the risk model

Using R software (“survivor”, “survminer” package), survival curves were plotted for patients of different genders, tumor grades, and tumor stages based on the risk model to verify the applicability of our constructed risk model to patients of different clinical trait groups.

Functional and pathway enrichment analysis

GO function enrichment analysis: Using R software (“colorspace”, “stringi”, “ggplot2”, “circlize”, “RColorBrewer”, and “ggpubr” packages), risk-differentiated Disulfideptosis-associated lncRNAs were analyzed for which functions they were enriched with and visualized, including cellular component (CC), molecular function (MF), and biological process (BP), presented as bar and circle charts. KEGG pathway enrichment analysis: Using R software (“colorspace”, “stringi”, “ggplot2”, “circlize”, and “RColorBrewer” packages), KEGG pathway enrichment analysis was performed for risk-differentiated Disulfideptosis-associated lncRNAs and visualized as bubble graphs and circle plots.

Analysis of tumor-infiltrating immunocyte and immune checkpoints

GSVA analysis was performed using R software (“limma”, “GSVA”, “GSEABase”, “pheatmap”, and “reshape2” packages) on samples from high and low-risk groups, looking for correlations with immune function, and presented as heat maps.

Tumor immune dysfunction and exclusion (TIDE) and estimation of tumor mutational burden (TMB)

Based on the expression profile of TCGA, TIDE data of all samples were downloaded from TIDE (http://tide.dfci.harvard.edu/), and TIDE values were calculated using R software (“limma”, “ggpubr” package) to calculate the TIDE values of samples from high and low-risk groups and compare them to assess the immune escape ability of different high and low-risk groups and visualize them as box plots.

Using the R software (“maftools” package), the TMBs in the high and low-risk groups were analyzed based on the simple nucleotide variation data downloaded from TCGA above, the differences were compared, and the different mutation types were displayed, and the 15 genes with the highest mutation frequencies were selected for display and visualized in waterfall plots. The relationship between different high and low TMB and patient survival time was analyzed using R software (“survival”, “survminer” package), and KM curves were generated. Also, the relationship between different high and low TMB in the high and low-risk groups and patient survival time was analyzed and KM curves were generated.

RNA extraction and RT-qPCR

The RNA from cell lines and tissues was extracted using Trizol (Shanghai Kangwei Biotech, China), chloroform, isopropanol, 75% ethanol (Solarbio, Beijing, China), and RNase-free water. The extracted RNA was quantified, followed by reverse transcription into cDNA using a reverse transcription kit (Servicebio, Wuhan, China). Subsequently, real-time fluorescence quantitative PCR was performed using 2 × SYBR Green PCR Master Mix (Servicebio, Wuhan, China). The relative expression levels of the target genes were calculated using the 2^−ΔΔCt method. The experiments were repeated three times using different PCR runs, and the average values were obtained. β-actin was used as the internal reference gene.

Primer sequences

β-actin F: CATGTACGTTGCTATCCAGGC R: CTCCTTAATGTCACGCACGAT. Human-FAM225B F: TCGGGATAGTGATGGCAAGC R: GCCACCTCTGTGGTCCTAAC. Human-LINC01338 F: GCCCAGGCTTCCCGATTATT R: ATGGGTTTGACCGTCGATGT. Human-SPINT1-AS1 F: CCGGGTACTTGAGCTCCCTA R: TGATCAGCCCGGGAGACTTT. Human-WWC2-AS2 F: GGGTCGTGTTTGCCCTTAGA R: GCCCTAAATGCGGTCAAAGC. Human-ZNF503-AS1 F: CATTCTCCACCCTGCCACAT R: AGGCATCTTGGCAGAAGGAC.

Cell transfection and cell function experiments

FAM225B has been studied in other tumors or diseases and is associated with the development of nasopharyngeal carcinoma, glioblastoma, etc. [12–15]. However, relatively few studies have been conducted in renal clear cell carcinoma, so FAM225B was selected as the target of the follow-up study, and the sequence design for the FAM225B siRNA was as follows:

• Si1: ATGCTTTCTGCAAAGAATAATAC

• Si2: TGGTGCATGCATCTTTTCTGC

• Si3: ATGTGTAATGAACATTTAAAATT

Using RNAi technology, transient transfection of siRNA was performed in renal cancer cell lines OS-RC-2 and 786-O using Lipofectamine 2000 (Invitrogen, USA), opti-MEM medium (Life Technologies, USA), and other reagents. After culturing the cells under optimal conditions, RNA was extracted from the cells, and RT-qPCR was performed to determine the knockdown efficiency. The above constructed cells were inoculated in 96-well plates, and the fluorescence values at 450 nm (OD value) of these two cells at 0 h, 24 h, 48 h, 72 h, and 96 h were measured using CCK-8 kit (Solarbio, Beijing, China), to compare the cell proliferation ability.

The logarithmic growth phase of the constructed transfected siRNA kidney cancer cells OS-RC-2 and 786-O were inoculated in six-well plates, and when the cells grew to 80–90% using a 1 ml pipette tip scratch, PBS (Phosphate Buffered Saline) was washed followed by the addition of serum-free culture medium, photographed by microscopy and recorded for 0 h for each subgroup, after which they were placed in a 37°C, 5% CO₂ incubator. The 6-well plates were subsequently removed at 24 h to aspirate and discard the old medium, and PBS was gently blown to wash away the floating cells before adding the serum-free medium and continuing to take photographs to record the photographs of each subgroup.

Liquid Matrigel matrix gel (Corning, USA) was diluted 1:9 using RIPM-1640 medium without FBS (Fetal Bovine Serum) and double antibodies (Anti-penicillin and streptomycin). Afterward, 70 ul of diluted substrate gel was taken and spread flat in Transwell chambers in 24-well plates and left in the incubator for 1 hour to allow the substrate gel to solidify. After digestion of the cells, the cells were resuspended using 200 ul of RIPM-1640 medium without FBS according to the cell density, and the cell volume was adjusted to 5 × 10⁴ cells and placed in the upper chamber. After adding 600 ul of RIPM-1640 complete medium in the lower chamber, they were placed in the cell culture chamber. The culture was terminated at 48 h, the chambers were removed, and the medium was aspirated and washed twice with PBS. The Matrigel matrix gel and cells in the chambers were wiped with cotton swabs, and 600 ul of methanol was added to the lower chamber for 30 min. 600 ul of Giemsa stain was added to the lower chamber after aspiration of the methanol, and the cells were stained for 30 min.

Statistical analysis

All data were statistically analyzed using R4.1.0 software (https://www.r-project.org/) and Kaplan-Meier survival analysis was used to analyze the differences between the two risk groups. Independent samples t-test was used to compare the differences between the groups and statistical difference was considered at p < 0.05.

Data availability statement

The data supporting this study’s findings are available from the corresponding author upon reasonable request.

Acquisition of differentially expressed disulfidptosis-related lncRNAs

From previous studies, we obtained 13 genes associated with Disulfideptosis. Through co-expression analysis, we identified 892 lncRNAs that showed co-expression relationships with these 13 Disulfideptosis-associated genes. The co-expression relationships are depicted in a Sankey diagram shown in Figure 1. Subsequently, we performed univariate Cox regression analysis on the aforementioned lncRNAs and identified 48 lncRNAs that were associated with the prognosis of renal clear cell carcinoma. The corresponding forest plot is shown in Figure 2A. After Lasso regression analysis and multifactorial Cox regression analysis (Figure 2B, 2C), 5 lncRNAs (FAM225B, ZNF503-AS1, SPINT1-AS1, WWC2-AS2, LINC01338) emerged as independent risk factors for the prognosis of renal clear cell carcinoma. The expression relationships of these five LncRNAs obtained with 13 genes related to Disulfidptosis were shown in Figure 2D.

Figure 1. Co-expression of disulfidptosis-associated genes and lncRNAs.

Figure 2. Screening of lncRNA for independent prognosis. One-way Cox regression analysis (A); Lasso regression analysis (B); λ values (C); correlation heat map (D).

Development of a prognostic risk model

According to the formula Risk Score = (FAM225B × (1.21052350282848)) + (ZNF503-AS1 × (−0.413072607271419)) + (SPINT1-AS1 × (−0.348131526553965)) + (WWC2-AS2 × (−0.289726017312162)) + (LINC01338 × (1.00289905657691)). This indicates that FAM225B and LINC01338 are highly expressed in the high-risk group, while ZNF503-AS1, SPINT1-AS1, and WWC2-AS2 are lowly expressed in the high-risk group. We randomly divided the 533 samples in TCGA into two groups (train group with 267 cases and test group with 266 cases). The t-test results for the clinical data of both groups (Age, Gender, Grade, Stage, T, M, N) are presented in Supplementary Table 1. From the table, it is evident that there were no significant differences in all clinical data between the two randomly assigned groups (all p > 0.05), indicating satisfactory randomness of the two groups.

Prognostic features of risk models

Based on the risk scores obtained from the previous risk model, patients were divided into high and low-risk groups using the median method. The risk curves, survival status curves, and risk heatmaps of all patients are shown in Supplementary Figure 1A–1C. It is evident that our risk model effectively stratifies patients into high and low-risk groups. As the risk score increases, the number of patient deaths also increases, which aligns with our expectations. From the risk heatmap, it is also evident that FAM225B and LINC01338 are positively correlated with the risk score, while ZNF503-AS1, SPINT1-AS1, and WWC2-AS2 are negatively correlated with the risk score. This pattern is consistent in both the randomly assigned train and test groups, demonstrating the consistent and significant performance of the risk model. This internal validation of the risk model is shown in Figure 3A–3F.

Figure 3. Risk curve. Train set (A–C); Test set (D–F).

Clinical features and risk scores

After dividing all samples into high and low-risk groups based on the risk scores, Kaplan-Meier survival curves were used to compare the overall survival rates and progression-free survival rates in the high and low-risk groups of all samples. It is evident that both the overall survival time (Figure 4A) and progression-free survival time (Figure 4B) of the high-risk group are significantly lower than those of the low-risk group (all p < 0.001). To further validate the findings, in both the train and test sets, the high-risk group exhibited significantly lower overall survival time (Figure 4C, 4D) compared to the low-risk group (all p < 0.001). This indicates that the risk model can effectively distinguish patients with different survival outcomes, demonstrating its favorable predictive performance.

Figure 4. Survival curve. OS (A) and PFS (B) of the total sample. OS of the train set (C) and the test set (D).

Nomogram and clinical indicators

To validate whether our constructed prognostic model is independent of other clinical characteristics (Age, Gender, Grade, Stage) as prognostic factors, we performed independent prognostic analysis. The results of univariate Cox regression analysis and multivariate Cox regression analysis showed that the risk score had a significant impact on prognosis, with p < 0.001 and hazard ratio of 1.079 (1.059–1.098) in univariate analysis, and p < 0.001 and hazard ratio of 1.043 (1.021-1.066) in multivariate analysis. These results are presented in Figure 5A, 5B, indicating that our constructed prognostic model is independent of other clinical characteristics and serves as an independent prognostic factor.

Figure 5. Independent prognostic analysis. Univariate analysis (A). Multivariate analysis (B). The ROC curve of the risk model for 1, 3, and 5 years (C). The 5-year ROC curve of the risk model in relation to other clinical characteristics (D).

Subsequently, we calculated the AUC (Area Under the Curve) values of the prognostic model and plotted the ROC (Receiver Operating Characteristic) curve. The AUC values for 1, 3, and 5-year survival were 0.756, 0.752, and 0.781, respectively, indicating a high predictive value (Figure 5C). The AUC value (Figure 5D) for the 5-year survival (0.781) is higher compared to other clinical characteristics such as Age (0.598), Gender (0.488), Grade (0.680), and Stage (0.717). This indicates that our constructed prognostic model has a significant advantage over individual clinical characteristics in terms of predictive accuracy.

Through Principal Component Analysis (PCA), we can observe whether the lncRNAs involved in model construction can differentiate patients in the high and low-risk groups. From the three-dimensional visualization of the PCA results (Figure 6A–6D), it is evident that in Figure 6D, patients in the high-risk group are positioned in the upper left region of the three-dimensional plot, while patients in the low-risk group are positioned in the lower right region, with a clear boundary. Therefore, the lncRNAs involved in model construction can effectively distinguish patients in the high and low-risk groups.

Figure 6. PCA analysis. All genes (A); Disulfidptosis genes (B); Disulfidptosis lncRNAs (C); Risk lncRNAs (D).

Based on the risk model and related clinical data, we constructed a nomogram (Figure 7A). We also generated calibration curves for the nomogram at 1, 3, and 5-year intervals (Figure 7B). It can be observed that the three curves are relatively close to the reference line, indicating that the nomogram can accurately predict the survival period of patients.

Figure 7. Construction and validation of nomogram column line plots constructed from risk models. A nomogram integrates risk scores and clinical characteristics to predict overall survival at 1, 3, and 5 years (A). Calibration curves for overall survival at 1, 3, and 5 years (B). ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

To investigate whether the risk model can be applied to patients in different clinical subgroups, we selected gender, Grade stage, and Stage classification as the clinical subgroups. From Figure 8A–8F, it is evident that our constructed risk model exhibits significant differences in survival probabilities between the high and low-risk groups in the FEMALE and MALE subgroups (Figure 8A, 8B), G1&2 and G3&4 subgroups (Figure 8C, 8D), and Stage I–II and Stage III–IV subgroups (Figure 8E, 8F) (all p < 0.001). This indicates that the risk model demonstrates excellent predictive performance in patients from different clinical subgroups as well.

Figure 8. Clinical subgroup analysis of the risk model. Gender (A, B). Tumor grading (C, D). Tumor staging (E, F).

Functional and pathway enrichment analysis in the risk model

The bar plot of the GO analysis (Figure 9A) reveals that in terms of Molecular Function (MF), the main focuses are on antigen binding and secondary active transmembrane transponder activity. In Cellular Component (CC), the main focuses are on the apical part of the cell and the collagen-containing extracellular matrix. In Biological Process (BP), the main focuses are on the defense response to bacterium and organic anion transport. GO enrichment analysis bubble plot (Figure 9B): The outermost circle represents the GO enrichment IDs, the second circle represents the number of genes associated with each GO term, the third circle represents the number of differentially expressed genes enriched in each GO term, and the fourth circle represents the proportion of enriched genes. The redder the color, the more significant the enrichment of differentially expressed genes. From the figure, we can observe that the gene count enriched in the GO term “0071735” is the highest in proportion. For KEGG pathway analysis, the bar plot (Figure 9C) indicates that the highest enrichment and differential expression are primarily observed in pathways such as Cytokine-cytokine receptor interaction, Neuroactive ligand-receptor interaction, and PI3K-Akt signaling pathway. From the circular plot of the KEGG pathway analysis (Figure 9D), it is also evident that the results are primarily focused on hsa00430 and hsa04610. From the GSEA analysis (Figure 9E), it can be observed that there are significant differences between the high and low-risk groups in terms of Type_II_IFN_Reponse, Type_I_IFN_Reponse, MHC_class_I, T_cell_co_stimulation, Inflammation-promoting, CCR, and Parainflammation. The calculated TIDE values (Figure 9F) are significantly lower in the low-risk group compared to the high-risk group (p < 0.0001). This finding further indicates that our predictive model can effectively distinguish the immune escape capabilities of patients in the high and low-risk groups.

Figure 9. Functional enrichment analysis. GO analysis (A, B). KEGG pathway analysis (C, D). GSEA analysis (E). Comparison of TIDE values between high and low-risk groups (F). ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Tumor mutational burden (TMB) analysis in the risk model

Based on the high and low TMB values, we divided the patients into high and low-risk groups (Figure 10A, 10B). The total mutation frequency in the low-risk group was 80.2%, while in the high-risk group, it was 78.92%, which is lower than that in the low-risk group. In the low-risk group, Nonsense Mutation was predominantly observed, while in the high-risk group, Missense Mutation was the major type. The top 15 genes with mutations, ranked by mutation frequency from high to low, are VHL, PBRM1, TTN, SETD2, BAP1, MTOR, MUC16, DNAH9, KDM5C, DST, LRP2, HMCN1, CSMD3, KMT2C, and FBN2. The KM curve of TMB in the high and low-risk groups shows that patients with low TMB have a higher probability of long-term survival. The combined analysis of TMB and the risk model score also reveals clear trends in the KM curve (Figure 10C). Patients with low TMB and low-risk scores tend to have better survival outcomes, while those with high TMB and high-risk scores have the poorest prognosis (Figure 10D).

Figure 10. Tumor mutational burden (TMB) analysis. Waterfall plot of TMB in high and low-risk groups (A, B). Overall survival analysis of TMB (C). Survival analysis of TMB in different high and low-risk groups (D).

Results of in vitro cell function experiments

The validated expression levels of the five selected lncRNAs associated with Disulfideptosis (FAM225B, ZNF503-AS1, SPINT1-AS1, WWC2-AS2, LINC01338) in five cell lines (HK-2, 786-O, OS-RC-2, ACHN, A-498) are shown in Figure 11A–11E. From the results, it is evident that FAM225B and WWC2- AS2 are highly expressed in multiple tumor cell lines. The interference efficiency of the two renal tumor cell lines (OS-RC-2, 786-O) is depicted in Figure 11F, 11G. Subsequently, Si1 and Si2 were selected for further studies.

Figure 11. Validation of 5 Disulfidptosis-associated lncRNAs. Comparison of RNA expression in renal cancer cells (786-O, OS-RC-2, ACHN, A-498) and normal renal proximal tubular epithelial cells (HK-2) (A–E). Validation of knockdown potency of FAM225B (F, G). Validation of cellular value added after FAM225B knockdown (H, I). ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, ^nsp > 0.05.

The OD values at 450 nm of OS-RC-2 and 786-O cells, measured by CCK-8 assay, after interference with Si1 and Si2 are shown in Figure 11H, 11I. It is evident that the proliferation of both cell lines is significantly reduced after interference.

In the in vitro cell migration assay, the wound healing capacity of OS-RC-2 cells after interference with Si1 and Si2 (Si1 and Si2 groups) is significantly lower compared to the control group (NC group) (Figure 12A, 12B). Similarly, this is also observed in 786-O cells (Figure 12C, 12D). These results indicate that the decreased expression of the LncRNA associated with Disulfideptosis (FAM225B) significantly affects the migration ability of both renal cancer cell lines.

Figure 12. Cellular function validation of FAM225B in two cell lines, 786-O, and OS-RC-2. Wound healing assay (A–D). Cell invasion assay (E–G). ^***p < 0.001.

In the cell invasion assay, the invasive capacity of OS-RC-2 cells after interference with Si1 and Si2 (Si1 and Si2 groups) is significantly lower compared to the control group (NC group) (Figure 12E, 12F). Similarly, this is also observed in 786-O cells (Figure 12G). These results indicate that the decreased expression of the lncRNA associated with Disulfidptosis (FAM225B) significantly affects the invasion ability of both renal cancer cell lines.