A Novel Ferroptosis-Related lncRNA Signature for Prognosis Prediction in Patients with Papillary Renal Cell Carcinoma

Introduction

Renal cell carcinoma (RCC) is one of the common malignancies originating from the renal tubular epithelium, with an incidence of approximately 431,288 in 2020, accounting for 2.2% of common malignancies worldwide.¹ Renal papillary cell carcinoma (PRCC) is the second most common renal cell carcinoma after renal clear cell carcinoma, accounting for 15%–20% of RCC cases.² On the basis of different clinicopathological and immunohistochemical features, PRCC can be divided into two subtypes: type 1 and type 2; in general, type 2 PRCC has a worse prognosis than type 1 PRCC.³ Although the comprehensive treatment of PRCC has improved, the overall survival (OS) rate remains low, especially for patients with advanced or metastatic tumors with limited treatment options.⁴ Therefore, new treatment targets and prognostic biomarkers for PRCC are necessary.

Unlike cellular autophagy and apoptosis, ferroptosis is an oxidative, iron-dependent form of necrotic cell death that is driven by the lethal accumulation of lipid peroxidation. Cancer cells have an increased iron requirement compared with normal non-cancerous cells, and studies have found that increased iron-dependent ROS-induced extra-mitochondrial lipid peroxidation triggers ferroptosis to inhibit tumorigenesis, making ferroptosis inducers a potential therapeutic agent for tumors.^5,6 Ferroptosis-related genes can be divided into three groups: driver, suppressor, and marker genes.⁷ Long non-coding RNAs (lncRNAs) are a subset of non-coding RNA molecules of about 200 nt, and studies have shown that lncRNAs play important roles in many life processes such as epigenetic regulation, cell cycle regulation, and cell differentiation regulation.^8–10 LncRNAs are involved in tumorigenesis, development, and metastasis processes.^11,12 Currently, studies on iron death-related lncRNAs are limited. Lu et al reported that lncRNA PVT1 regulates ferroptosis through miR-214-mediated TFR1 and p53.¹³ LncRNA LINC00336 suppresses ferroptosis in lung cancer by acting as a competitive endogenous RNA.¹⁴ The relationship between ferroptosis-related lncRNA and clinicopathological features of PRCC is unclear.

The Cancer Genome Atlas (TCGA) is one of the most ambitious and successful cancer genomics programs to date. The TCGA program has generated, analyzed, and provided available genomic sequence, expression, methylation, and copy number variation data on over 11,000 individuals who represent over 30 different types of cancer.¹⁵ In this study, we constructed the prognostic signature of differentially expressed ferroptosis-related lncRNAs in PRCC based on TCGA data. We confirmed that this lncRNA signature could be used as an effective and independent prognostic biomarker in patients with PRCC.

Materials and Methods

Data Collection and Preprocessing

The RNA-seq transcriptome data and clinical data of patients with PRCC were downloaded from TCGA database (https://portal.gdc.cancer.gov/). A total of 32 adjacent normal tissues and 289 cancer tissues from patients with PRCC were enrolled. The ferroptosis-related genes were downloaded from FerrDb (http://www.zhounan.org/ferrdb/)¹⁶, and we identified 259 ferroptosis-related genes as shown in Table 1.

Table 1 Ferroptosis-Related Genes from Downloaded from FerrDb

Identification and Differential Expression Analysis of Ferroptosis-Related lncRNAs

Pearson correlation analysis between the ferroptosis-related genes and lncRNAs was used to identify the ferroptosis-related lncRNAs (| R| > 0.5 and p < 0.001). The Wilcoxon signed-rank test was utilized to identify differentially expressed ferroptosis-related genes and differentially expressed ferroptosis-related lncRNAs by the limma package of R based on the cutoff values: |Log2FC| > 1 and FDR < 0.05.

The “clusterProfiler” package was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses on differentially expressed ferroptosis-related genes.^17,18

Development of the Ferroptosis-Related lncRNA Prognostic Signature

Lasso penalized Cox regression and univariate Cox regression analysis were employed to construct the ferroptosis-related lncRNA prognostic signature, which was used to construct a risk scoring function. The following formula was used to calculate the riskScore with the constructed risk model for all clinical cases: RiskScore = (coefficient lncRNA1 × expression of lncRNA1) + (coefficient lncRNA2 × expression of lncRNA2) + ⋯ + (coefficient lncRNAn × expression lncRNAn). The patients were classified in either low-risk (<median number) or high-risk (≥ median number) group based on the median score.

The Kaplan–Meier (K-M) survival curve described the predictive power of indicators using the “survival” and “survminer” packages. The area under the curve (AUC) of the receiver operating characteristic (ROC) curve and C-index were applied to verify the predictive accuracy using the “timeROC” and “rms” packages.

Building a Predictive Nomogram

On the basis of all independently predictive variables, we constructed a prediction nomogram using the “rms” and “foreign” packages for the prediction of 1, 3, and 5-year OS of patients with PRCC.

Identification of the Relationship Between the Risk Score and the Clinicopathological Characteristics

To better understand the role of risk score in the clinical outcome of PRCC, we systematically investigated the relationship between risk score and clinicopathological characteristics (including age, sex, stage status, T status, M status, and N status), as well as the expression level of prognostic lncRNAs, using the “pheatmap” package.

Protein–Protein Interaction Network Construction

The String database (https://cn.string-db.org/) was used to analyze the interactions of ferroptosis-related genes and ferroptosis-related lncRNAs. Cytoscape¹⁹ was used for the visualization of the protein–protein interaction (PPI) network.

Gene Set Enrichment Analyses

Gene set enrichment analysis (GSEA)²⁰ was performed to explore the potential altered pathways between the high- and low-risk groups.

Immunity Analysis

The CIBERSORT,^20,21 ESTIMATE,²² MCPcounter,²³ single-sample gene set enrichment analysis (ssGSEA),²⁴ and TIMER²⁵ algorithms were used to assess cellular components or cell immune responses. A total of 22 immune cell types were used to estimate relative abundance of immune cell infiltration between the low- and high-risk groups. ssGSEA was then used to quantify tumor-infiltrating immune cell subgroups during immune responses and immune functions between the two groups. We also compared the expression levels of immune checkpoints between the two groups.

Chemotherapy Sensitivity Analysis

We obtained the drug sensitivity data and RNA-seq data of the NCI-60 cell line from CellMiner (https://discover.nci.nih.gov/cellminer/home.do) and analyzed the relationship between the prognostic lncRNAs and drug sensitivity using Pearson correlation analysis.

Statistical Analysis

R version 4.0.2 software and its resource packages were used for statistical analysis and related visualization graphics. Student’s t test was used to compare differences between different risk groups. Statistical significance was set at p <0.05.

Results

Enrichment Analysis of Ferroptosis-Related Genes

The work flowchart is shown in Figure 1. The detailed clinical characteristics of these patients are shown in Table 2. We identified 56 differentially expressed ferroptosis-related genes (32 upregulated and 24 downregulated, Table S1). Differentially expressed genes (DEGs) were enriched in the biological processes with response to oxidative stress, cellular response to chemical stress, and cellular response to oxidative stress. They were enriched in molecular functions with oxidoreductase activity, acting on NAD(P)H, oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, and iron ion binding. Cellular component (CC) was mainly enriched in apical parts of the cell, oxidoreductase complex, and NADPH oxidase complex. The biological pathways in KEGG were enriched in arachidonic acid metabolism, biosynthesis of amino acids, leishmaniasis, HIF-1 signaling pathway, and carbon metabolism (Figure 2).

Table 2 The Clinical Characteristics of Patients with PRCC in the TCGA Dataset

Figure 1 Flowchart of this study.

Figure 2 GO and KEGG analyses of differentially expressed ferroptosis-related genes. (A) Functional annotation of differentially expressed ferroptosis-related genes between adjacent normal and tumor samples using GO terms. (B) Functional annotation of differentially expressed ferroptosis-related genes between adjacent normal and tumor samples using KEGG terms.

Ferroptosis-Related lncRNAs’ Prognosis Signature in PRCC

A total of 680 ferroptosis-related lncRNAs were identified (Table S2). Univariate Cox analysis was performed to evaluated the relationship between the expression of ferroptosis-related lncRNAs and patients’ OS, and the results showed that 29 ferroptosis-related lncRNAs (DARS-AS1, AL355512.1, ZFAS1, LUCAT1, AL022344.1, LINC01559, AC010624.2, AL031710.1, AL355102.4, MNX1-AS1, AC109460.1, AC127537.1, CD44-AS1, AC099850.4, LINC02154, AC024022.1, AC092115.4, AC026401.3, CYTOR, AL365181.3, LINC02535, ADAMTS9-AS1, AC022075.1, AC096642.1, AC107464.2, SLC25A5-AS1, AC104564.1, ITGB1-DT, and MIR4435-2HG) were independent prognostic predictors of PRCC (Figure 3). Multivariate Cox analysis indicated that 15 of 29 ferroptosis-related lncRNAs were the best candidates to construct the prognostic signature (Table 3). The risk score was calculated for each patient with PRCC as follows: −0.0152169179410617 × the expression levels of ZFAS1 + (−0.814201526659232) × the expression levels of AC010624.2 + (−0.574312213529718) × the expression levels of AL031710.1 + (−0.295711007986573) × the expression levels of AL355102.4 + 0.343550090287168 × the expression levels of MNX1-AS1 + 2.59391370444094 × the expression levels of AC109460.1 + 0.514983257942382 × the expression levels of AC127537.1 + 0.178442445321917 × the expression levels of AC099850.4 + 0.0362308172908584 × the expression levels of LINC02154 + (−3.2728034303393) × the expression levels of AC024022.1 + 0.142403529727197 × the expression levels of AC026401.3 + 0.0664305502643028 × the expression levels of LINC02535 + (−0.0190675805173939) × the expression levels of ADAMTS9-AS1 + (−0.194060378575192) × the expression levels of AC107464.2 + 0.367056382643099 × the expression levels of MIR4435-2HG. Thus, we constructed a ferroptosis-related lncRNA prognosis model to predict the prognosis of patients with PRCC.

Table 3 Fifteen Prognostic Ferroptosis-Related lncRNAs Identified by Multivariate Cox Analysis

Figure 3 Univariate Cox regression analysis identified 29 ferroptosis-related lncRNAs as independent prognosis predictors of PRCC.

Survival Results and Multivariate Examination

All patients were divided into the high- and low-risk groups based on the median values of risk scores. K-M analyses revealed that the OS of the high-risk group was significantly worse than that of the low-risk group (p<0.001, Figure 4A). The AUC was 0.930, indicating its excellent predictive performance (Figure 4B). The AUC scores at 1, 2, and 3 years were 0.930, 0.953, and 0.933, respectively (Figure 4C). The survival status plot showed that the patients’ risk score was inversely proportional to the survival of patients with PRCC (Figure 4D and E). DCA also indicating the excellent predictive performance of the risk model (Figure 4F). The heatmap demonstrated that 10 lncRNAs were positively correlated with the risk model (Figure 4G).

Figure 4 Ferroptosis-related lncRNAs’ prognostic signature in patients with PRCC. (A) Kaplan–Meier survival curve result. (B) AUC values of the risk factors. (C) AUC values at 1, 2, and 3 years (D) Distribution of risk scores of patients with high- and low-risk PRCC based on the ferroptosis-related lncRNAs’ prognostic signature. (E) Scatter plot shows the correlation between the survival status of the high- and low-risk groups. (F) DCA of the risk factors. (G) Heat map of the expression of 15 ferroptosis-related lncRNAs between the high- and low-risk groups.

Univariate and multivariate Cox analyses were used to test whether risk score is an independent prognostic factor for OS. The lncRNA signature (HR:1.005, 95CI: 1.002–1.007) and tumor stage (HR:2.232, 95CI: 1.640–3.038) were independent prognostic factors of OS of patients with PRCC (Figure 5A and B). The hybrid nomogram incorporating clinicopathological characteristics and the novel ferroptosis-related lncRNAs’ prognostic signature might be applied in the clinical management of patients with PRCC (Figure 5C).

Figure 5 Univariate (A) and multivariate COX (B) analyses between overall survival rate and various clinicopathological parameters. (C) Nomogram for both clinicopathological factors and prognostic ferroptosis-related lncRNAs (*p <0.05; **p <0.01; ***p <0.001).

The association between ferroptosis-related lncRNAs’ prognostic signature and clinicopathological manifestations was shown in the heatmap (Figure 6).

Figure 6 Heatmap for ferroptosis-related lncRNAs’ prognostic signature and clinicopathological manifestations (*p <0.05; **p <0.01; ***p <0.001).

Gene Set Enrichment Analysis

The relationship between ferroptosis-related genes and lncRNAs is shown in Figure 7. GSEA revealed that the prognostic ferroptosis-related lncRNAs were associated with the tumor-related pathways and immune-related pathways, such as acute myeloid leukemia, bladder cancer, chemokine signaling pathway, primary immunodeficiency, natural killer cell-mediated cytotoxicity, and T cell receptor signaling pathway (Figure 8).

Figure 7 Relationship between the prognostic ferroptosis-related lncRNA and mRNA expression levels.

Figure 8 Gene enrichment analysis for prognostic ferroptosis-related lncRNAs.

Immune Status Difference Between Patients with High-Risk and Low-Risk PRCC

The immune response heatmap based on CIBERSORT, ESTIMATE, MCP counter, ssGSEA, and TIMER algorithm is shown in Figure 9A. ssGSEA was used to analyze the correlation between immune cell subsets and related functions, and the results showed significant differences in inflammation promotion, parainflammation, and Type I IFN response between the low- and high-risk groups (p <0.01; Figure 9B). In addition, the expression of immune checkpoints including CD80, IDO1, and LAG3 was significantly different between the low-risk and high-risk groups (Figure 9C).

Figure 9 (A) Immune status difference between high-risk and low-risk groups. Heatmap for immune responses based on CIBERSORT, ESTIMATE, MCPcounter, ssGSEA, and TIMER algorithms among high- and low-risk groups. (B) Immune-related function between high- and low-risk groups. (C) Expression of immune checkpoints between high- and low-risk groups. *p <0.05; **p <0.01; ***p <0.001.

Prognostic lncRNA Expression and Chemotherapy Sensitivity

We investigated the expression of prognostic lncRNAs in NCI-60 cell lines and analyzed the relationship between the expression levels and drug sensitivity. The results showed that MNX1-AS1, ZFAS1, MIR4435-2HG, and ADAMTS9-AS1 were significantly correlated with the sensitivity of some chemotherapy drugs (Figure 10 and Table 4). For instance, MNX1-AS1 could increase the sensitivity of most chemotherapy drugs, including fulvestrant, LDK-378, fluorouracil, tegafur, and PF-06463922. ZFAS1 could reduce the sensitivity of most chemotherapy drugs, including everolimus and rapamycin.

Table 4 Prognostic lncRNAs Expression and Chemotherapy Sensitivity

Figure 10 Prognostic lncRNA expression and chemotherapy sensitivity.

Discussion

Given that the molecular mechanisms of PRCC and clear cell renal cell carcinoma (RCCC) differ, and the incidence of PRCC is lower than that of PCCC, patients with PRCC are always excluded from large clinical trials of targeted drugs, and research on PRCC is not as popular as studies on RCCC.²⁶ Although some patients with PRCC can receive early surgical treatment, a large number of advanced patients continue to die due to postoperative recurrence, tumor metastasis, or resistance to chemotherapy drugs. Therefore, new and effective prognostic biomarkers for PRCC should be explored.²⁷ This study is the first to establish a prognostic model based on ferroptosis-related lncRNAs in patients with PRCC.

In this study, we first identified 680 ferroptosis-related lncRNAs, analyzed 29 prognostic ferroptosis-related lncRNAs via multivariate Cox regression analysis, and established a ferroptosis-related lncRNA signature model for PRCC. The prognostic model of PRCC proposed in this study consisted of 15 ferroptosis-related lncRNAs. Studies have reported that these lncRNAs play important roles in the pathogenesis of tumors. ZFAS1 is a major regulator of the epithelial–mesenchymal transition (EMT) through miR-200/ZEB1/E-cadherin and vimentin signaling in colon adenocarcinoma.²⁸ The expression level of AL031710.1 decreases with the increase in tumor stage of papillary renal cell carcinoma.²⁹ AL355102.4 is downregulated in clear cell renal cell carcinoma samples and indicates good prognosis when lowly expressed, which reveals that AL355102.4 is a risk factor with high gene expression.³⁰ High MNX1-AS1 levels are correlated with poor survival outcome and aggressive clinicopathological characteristics in various cancers, suggesting that MNX1-AS1 may be applied as a prognostic marker and even a therapeutic target.³¹ In the miRNA–lncRNA pair centered regulatory network consisting of competing triplets in ovarian cancer, AC099850.4 is one of the top five lncRNAs with the most connections.³² Studies have shown that LINC02154 is significantly upregulated in renal cell carcinoma, and its high expression is one of the risk factors for poor prognosis. It is involved in the construction of a model for predicting the prognosis of patients with renal cell carcinoma.^33,34 AC024022.1 is a protective biomarker in papillary renal cell carcinoma.³⁵ LINC02535 co-functions with PCBP2 to regulate DNA damage repair in cervical cancer by stabilizing RRM1 mRNA.³⁶ LINC02535 is highly expressed in lung Adenocarcinoma tissues and closely related to the OS of lung adenocarcinoma. Knocking down LINC02535 also significantly inhibits the EMT, a key contributor to tumor invasion and metastasis.³⁷ ADAMTS9-AS1 restrains the aggressive traits of breast carcinoma cells via sponging miR-513a-5,³⁸ but it promotes cell proliferation and EMT in colorectal cancer.³⁹ MIR4435-2HG is involved in the genesis and progression of a variety of tumors.^40–42 Although these lncRNAs are closely related to the pathogenesis of tumors, we combined them as prognostic markers in patients with PRCC for the first time.

Our results demonstrated that this model had good predictive performance using Kaplan–Meier and ROC curves. The results of univariate and multivariate Cox analyses suggested that this lncRNA signature and tumor stage were independent prognostic factors of OS of PRCC. The hybrid nomogram incorporating clinicopathological characteristics and the novel ferroptosis-related lncRNAs’ prognostic signature might be applied in the clinical management of patients with PRCC.

As ferroptosis is associated with immunity, we evaluated immune cell infiltration and immune function between the low- and high-risk groups. Compared with the low-risk group, the high-risk group had a greater degree of immune cell infiltration, such as B cells, T cell CD8, macrophages, and NK cells, indicating that the high-risk group had stronger anti-tumor immune activity than the low-risk group. In addition, eeGSEA confirmed that the inflammation-related immune function in the high-risk group was more active than that in the low-risk group. We tested the expression level of immune checkpoints and found that the expression levels of most immune checkpoints were significantly higher in the high-risk group than in the low-risk group. Given that the expression of immune checkpoints is usually positively correlated with immunotherapy responsiveness, the high-risk group might better respond to immunotherapy than the low-risk group.

On the basis of the results of the chemotherapy drug sensitivity analysis, we found that the expression of some prognostic lncRNAs was related to the resistance or sensitivity of many Food and Drug Administration-approved chemotherapy drugs. The increased expression of MNX1-AS1 was associated with increased drug sensitivity. ADAMTS9-AS1 was negatively correlated with the sensitivity of Axitinib but positively correlated with sensitivity of Bosutinib, Erlotinib, Lapatinib, and Neratinib. Therefore, some prognostic lncRNAs can be used as therapeutic targets to overcome chemotherapy drug resistance or adjuvant chemotherapy drug sensitivity.

This study still had some limitations. We performed retrospective analysis using data from a public database, so bias was inevitable. Therefore, large prospective studies are needed to confirm the validity of the prognostic features.

Conclusion

We constructed and validated a ferroptosis-related lncRNA prognosis signature in PRCC, which consisted of 15 lncRNAs. This novel ferroptosis-related lncRNA prognosis signature could predict the survival of patients with PRCC. Besides, we found significant differences in immune cell infiltration and immune-related function between the high- and low-risk groups, which might help for follow-up immune intervention therapy. Thus, our study might provide a theoretical basis for promising biomarkers and targets for the individualized therapy of PRCC.