Abstract
Background
Ferroptosis is a novel type of programmed cell death in various tumors;
however, underlying mechanisms remain unclear. We aimed to develop
ferroptosis-related long non-coding RNA (FRlncRNA) risk scores to
predict lower-grade glioma (LGG) prognosis and to conduct functional
analyses to explore potential mechanisms.
Methods
LGG-related RNA sequencing data were extracted from The Cancer Genome
Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases. Pearson
correlation analysis was used to identify the FRlncRNAs, univariate Cox
regression analysis was for identify the prognostic FRlncRNAs, and then
intersection FRlncRNAs were screened between TCGA and CGGA. Least
absolute shrinkage and selection operator (LASSO) Cox regression was
used to develop a risk score to predict LGG prognosis.
Results
A total of nine FRlncRNAs were screened to construct the novel
prognostic risk score of LGG, and high-risk score patients had a worse
overall survival than low-risk score patients both in TCGA and CGGA
datasets. The risk score was quite correlated with clinicopathological
characteristics (age, WHO grade, status of MGMT Methtlation, IDH
mutation, 1p/19q codeletion, and TMB), and could promote current
molecular subtyping systems. Comprehensive analyses revealed that
signaling pathways of B-cell receptor and T-cell receptor, immune cells
of macrophage cell and CD4+ T cell, tumor microenvironment of stroma
score and immune score, and immune checkpoints of PD-1, PD-L1, and
CTLA4 were all enriched in the high-risk score group.
Conclusion
The nine FRlncRNAs risk scores was a promising biomarker to predict the
LGG’s prognosis and distinguish the characteristics of molecular and
immune.
Supplementary Information
The online version contains supplementary material available at
10.1007/s12672-024-01587-9.
Keywords: Ferroptosis, Long non-coding RNAs, Lower-grade glioma, Risk
score, Immune checkpoints
Introduction
Approximately 30,000 patients were diagnosed with glioma in the United
States from 2013 to 2017, accounting for approximately 81.1% of all
malignant central nervous system tumors [[32]1]. Surgical excision
followed by temozolomide-based adjuvant treatment is standard treatment
for gliomas, but long-term prognosis remains far from satisfactory
[[33]2]. The reasons for this poor prognosis may be (1) patients in the
same subtype stratified by the existing molecular subtyping often have
different outcomes [[34]3], (2) total removal of intracranial lesions
is anatomically difficult to perform [[35]2], and (3) while adjuvant
treatments based on the current molecular pathology have failed in
several trials [[36]4–[37]6]. Therefore, the investigation of novel
predictive biomarkers and therapeutic targets to treat gliomas is
urgently needed.
Ferroptosis is a newly discovered type of non-apoptotic programmed cell
death that is characterized by iron-dependent lipid peroxidation
[[38]7]. Previous findings suggest that ferroptosis can suppress tumor
growth and progression, and survival models based on
ferroptosis-related genes (FRGs) have been successfully constructed to
predict the prognosis of various types of malignant tumors
[[39]8–[40]10], including lower-grade glioma (LGG) [[41]11]. In
addition, ferroptosis is also associated with therapy resistance to
chemotherapy or radiotherapy, which provides insight into novel
therapeutic alternatives for glioma treatment, such as using iron
chelators [[42]12] and lipid peroxidation inhibitors [[43]13]. Although
FRGs are reported to be associated with the prognosis of gliomas due to
cell membrane destruction regulated by iron-dependent lipid
peroxidation, the transcription and translation of FRGs remain unclear
[[44]14].
Long non-coding RNAs (lncRNAs) are a subset of RNA molecules of
approximately 200 nt that play significant roles in gene transcription,
microRNA (miRNA) regulation, and protein folding [[45]15]. Preliminary
studies have found that aberrant lncRNA expression is associated with
glioma pathogenicity, such as that of lncRNAs HOXA11-AS [[46]16] and
MALAT1 [[47]17], and previous studies identified the stability of
lncRNA-based prognostic signatures in various kinds of cancers [[48]18,
[49]19]. Meanwhile, lncRNAs also participate in the regulation of
ferroptosis. Silencing of lncRNAs LINC00618 [[50]20], PVT1 [[51]21],
and ZFAS1 [[52]22] has been found to attenuate ferroptosis through
lipid peroxidation. However, there are currently no studies on
ferroptosis-related lncRNAs in gliomas.
In the current study, we developed a multi-ferroptosis-related lncRNAs
(FRlncRNAs) risk score to predict overall survival (OS) of LGG patients
using data from The Cancer Genome Atlas (TCGA), which was also
validated using the Chinese Glioma Genome Atlas (CGGA) databases. Then
comprehensive analysis was performed to explore the potential value of
the risk score, and a competing endogenous RNA (ceRNA) network was
developed to search for a potential target for LGGs.
Methods
Data sources
LGG in this study was defined as grade II or grade III gliomas
according to the World Health Organization (WHO) classification. RNA
sequencing data related to LGG in the TCGA database were extracted
using the Genomic Data Commons Data Portal
([53]https://portal.gdc.cancer.gov/). Related clinicopathological data
were acquired from the cBioPortal website
([54]https://www.cbioportal.org/). RNA expression files in the CGGA
database and their corresponding clinicopathological data were
downloaded from the CGGA website ([55]http://www.cgga.org.cn/). The
[56]GSE108474 and [57]GSE43378 datasets were retrieved from the Gene
Expression Omnibus (GEO) database ([58]www.ncbi.nlm.nih.gov/geo).
Clinical characteristics of the patients in the datasets are shown in
Table [59]1
Table 1.
Clinical characteristics of the LGG patients used in this study^a
TCGA cohort CGGA cohort [60]GSE108474 [61]GSE43378
No. of patients 515 625 187 18
Age
Median (range) 41 (14–87) 40 (10–74) 0 (0.0%) 51.5 (24–73)
≤ 40 (%) 254 (49.3%) 330 (52.8%) 0 (0.0%) 6 (33.3%)
> 40 (%) 261 (50.7%) 294 (47.0%) 0 (0.0%) 12 (66.7%)
Unknown 0 (0.0%) 1 (0.2%) 187 (100.0%) 0 (0.0%)
Gender (%)
Female 230 (44.7%) 263 (42.1%) 60 (32.1%) 4 (22.2%)
Male 285 (55.3%) 362 (57.9%) 94 (50.3%) 14 (77.8%)
Unknown 0 (0.0%) 0 (0.0%) 33 (17.6%) 0 (0.0%)
Grade (%)
WHO I 0 (0.0%) 0 (0.0%) 2 (1.1%) 0 (0.0%)
WHO II 249 (48.3%) 291 (46.6%) 92 (49.2%) 5 (27.8%)
WHO III 265 (51.5%) 334 (53.4%) 69 (36.9%) 13 (72.2%)
Unknown 1 (0.2%) 0 (0.0%) 24 (12.8%) 0 (0.0%)
MGMT methylation status (%)
Methylated 93 (18.1%) 298 (47.7%) NA NA
Unmethylated 422 (81.9%) 211 (33.8%) NA NA
Unknown 0 (0.0%) 116 (18.5%) NA NA
IDHMutant (%)
Mutant 414 (80.4%) 439 (70.2%) NA NA
Wild 92 (17.9%) 144 (23.1%) NA NA
Unknown 9 (1.7%) 42 (6.7%) NA NA
1p/19q Coldel (%)
Codel 173 (33.6%) 191 (30.6%) NA NA
Non-codel 210 (40.8%) 393 (62.9%) NA NA
Unknown 132 (25.6%) 41 (6.6%) NA NA
[62]Open in a new tab
a. Inclusion Criteria: 1. Gene expression profiles of patients with
grade I-III glioma
Exclusion Criteria: 1. Patients with grade IV glioma. 2. Duplicate
reports from the same cohort. 3. Protocol, case reports, or reviews. 4.
Data unavailable
lncRNAs annotation
The lncRNA annotation file (Genome Reference Consortium Human Build 38,
GRCh38) was downloaded from the GENCODE website
([63]https://www.gencodegenes.org/human/) and used to annotate the
lncRNA in the TCGA dataset. Accordingly, 14,247 lncRNAs were identified
in the TCGA dataset and 4304 identified in the CGGA dataset.
FRlncRNAs identification
A list of recognized FRGs, including 150 driver genes, 109 suppressor
genes, and 123 marker genes, was acquired from the FerrDb database
(Table S1) [[64]23]. Pearson correlation analysis (|R|> 0.3,
P < 0.0001) was used to identify FRlncRNAs in each dataset. The
FRlncRNAs in the two datasets were screened using univariate Cox
regression analysis based on OS (P < 0.001). The FRlncRNAs that
intersected both the TCGA and CGGA datasets were identified as
candidate FRlncRNAs for determining the risk scores.
FRlncRNAs risk score construction
Least absolute shrinkage and selection operator (Lasso)-penalized Cox
regression was used to construct the FRlncRNAs risk scores following
the minimum criteria with the penalty parameter (λ) based on “glmnet” R
package using a tenfold cross-validation (Fig. S1). The risk score of
each patient was calculated according to the following formula:
[MATH: Riskscore=∑i=1
nCoefi×xi :MATH]
Patients in each dataset were then divided into high-risk and low-risk
subgroups according to their risk scores. The median risk score was
used as the cut-off value. The risk score was evaluated using the
receiver operating characteristic (ROC) curve and validated using
univariate & multivariate Cox analyses based on OS in CGGA,
[65]GSE108474, and [66]GSE43378 datasets.
Function analysis
Gene ontology (GO) function analysis and Kyoto Encyclopedia of Genes
and Genomes (KEGG) pathway enrichment analysis were preformed according
to the differentially expressed genes (DEGs) between the high-risk and
low-risk groups. The DEGs were defined according to an absolute
log-twofold-change (|log2FC|) ≥ 1 and a false detection rate
(FDR) < 0.05. Gene set enrichment analysis (GSEA) was used to identify
potential biological processes and pathways that differed between the
high-risk and low-risk subgroups [[67]24]. The correlation between
estimated numbers of immune cells in the tumor environment and risk
scores was evaluated using various tools available on the TIMER 2.0
website ([68]http://timer.cistrome.org), including XCELL, TIMER,
QUANTISEQ, MCPcounter, EPIC, and CIBERSORT.
Hierarchical clustering
To examine the relationship between the newly identified molecular
subtypes and the previously characterized LGG immune subtypes, we
conducted unsupervised hierarchical clustering utilizing the CIBERSORT
algorithm, which quantifies the immune cell composition of each sample.
This analysis enabled us to categorize the patients into two distinct
subgroups, cluster 1 (n = 771) and cluster 2 (n = 369).
Analysis of tumor mutation burden (TMB)
Mutation data were also downloaded from the TCGA database, including
somatic mutations, insertion-deletion mutations, coding mutations, and
base replacement per million bases data. The TMB score was calculated
as the total number of somatic mutations. The LGG patients were then
separated into low-TMB and high-TMB groups using the median TMB value
as a cut-off value.
Predictive nomogram construction
Independent prognostic factors were identified using multivariate Cox
regression analysis. Based on the analysis results, a nomogram was
developed to predict the prognosis of patients with LGG.
Ferroptosis-related ceRNA network construction
To investigate the potential regulatory mechanism of FRlncRNAs in LGG
on mRNA expression via the sponging of miRNAs, the ceRNA network was
constructed. Briefly, FRlncRNA-related miRNAs were identified using the
miRcode database ([69]http://www.mircode.org/) and mRNAs related to the
miRNAs were identified using the miRDB ([70]http://mirdb.org/),
miRTarBase ([71]http://mirtarbase.mbc.nctu.edu.tw/php/index.php), and
TargetScan databases ([72]http://www.targetscan.org/). The ceRNA
network was plotted using Cytoscape 3.7.2 software [[73]25].
Statistical analysis
The Pearson correlation analysis used to evaluate the correlation
between lncRNAs and FRGs was performed using a |R^2|< 0.3 and P-value
of 0.0001. Kaplan–Meier analysis was used to compare OS between
different groups using the log-rank test. Univariate and multivariate
Cox regression analyses were performed to identify independent
predictors of OS. The diagnostic value of the risk scores and nomogram
model were evaluated using ROC curves. The Spearman test was used to
determine the association between estimated immune cell numbers and
risk scores. The Wilcoxon rank sum test was used to evaluate
differences stratified by various clinicopathological characteristics.
All statistical analyses were conducted using R software (version
4.0.3, [74]https://www.r-project.org). Unless specified otherwise, P
with two tailed value of less than 0.05 was considered statistically
significant.
Results
Identification of FRlncRNAs in LGG patients
Totally, 515 patients and 625 patients were screened in the TCGA
database, and CGGA database, respectively. A total of 1637 FRlncRNAs
were identified in the TCGA dataset and 493 in the CGGA dataset. In
addition, 667 prognostic lncRNAs in the TCGA dataset and 215 in the
CGGA dataset were confirmed to be associated with OS in LGG patients
(Table S2). Finally, 81 lncRNAs were identified as intersecting
prognostic FRlncRNAs between the two datasets (Table S3). A flow chart
of this study is depicted in Fig. [75]1a.
Fig. 1.
[76]Fig. 1
[77]Open in a new tab
a Flow chart of the study. b The competing endogenous RNA (ceRNA)
network of 31 ferroptosis-related lncRNAs (FRlncRNAs; red), 30 target
microRNAs (miRNAs; yellow), and 49 mRNAs (blue)
Ferroptosis-related ceRNA network construction
The ceRNA network is a complex regulatory network between non-coding
RNAs and coding RNAs in the cell via miRNAs. Of the 81 intersecting
FRlncRNAs, 31 were found in the miRcode database. Thirty of the miRNAs
were identified as being associated with FRlncRNAs with 280 pairs.
Subsequently, 49 mRNAs were found to be associated with 80 pairs of the
FRlncRNA-related miRNAs. The ceRNA network of the most critical 31
lncRNAs, 30 miRNAs, and 49 mRNAs is depicted in Fig. [78]1b.
Construction of a prognostic risk score in the TCGA training set
Using LASSO Cox regression analysis, nine FRlncRNAs, including
AGAP2-AS1, FAM181A-AS1, CRNDE, PAXIP1-AS2, TMEM254-AS1, SLC25A21-AS1,
DGCR9, LINC00844, and LINC00641, were screened to establish a risk
score for predicting the prognosis of LGG patients in the TCGA training
set, which is shown in Fig. [79]2a along with each coefficient value.
The prognostic risk score of each patient was determined according to
the following formula:
Fig. 2.
[80]Fig. 2
[81]Open in a new tab
a Coefficient values of nine ferroptosis-related lncRNAs (FRlncRNAs)
associated with a risk model for lower-grade glioma (LGG). b Sankey
diagram showing the top five mRNAs and risk type related to each
FRlncRNA. Kaplan–Meier curves (c), risk plot distribution and survival
status (d), and time-dependent receiver operating characteristic (ROC)
curves (e) of The Cancer Genome Atlas (TCGA) dataset. Kaplan–Meier
curves (f), risk plot distribution and survival status (g), and
time-dependent ROC curves (h) of the Chinese Glioma Genome Atlas (CGGA)
dataset
[MATH: Riskscore=AGAP2-<
/mo>AS1×0.066+FAM181A
-AS1×0.087+CRNDE×<
/mo>0.264+PAXIP1<
/mn>-AS2×0.330-TMEM254
-AS1×0.190-SLC25A2
1-AS1×0.109<
/mrow>-DGCR9×<
/mo>0.077-LINC00844×0.076-LINC00641×0.066 :MATH]
Nine FRlncRNAs related mRNAs were identified based on |R^2|< 0.3 and a
P-value of 0.0001 and the top five mRNAs related to each FRlncRNA are
shown in Fig. [82]2b. Of note, TMEM254-AS1, SLC25A21-AS1, DGCR9,
LINC00844, and LINC00641 were protective factors for OS, while
AGAP2-AS1, FAM181A-AS1, CRNDE, and PAXIP1-AS2 were risk factors for OS.
The median risk score of the patients in the TCGA dataset was − 0.754
(− 1.433–1.556). The patients were divided into low-risk and high-risk
groups based on the median risk score values and the groups exhibited
apparent differences in terms of OS (P < 0.001, Fig. [83]2c). The
corresponding risk scores and survival status are plotted in
Fig. [84]2d. The time-dependent ROC curves revealed the FRlncRNAs
exhibited better discrimination than the other prognostic models in the
TCGA dataset (age, gender, grade, and status of MGMT methtlation, IDH
mutation, 1p/19q codeletion; Fig. [85]2e). Similar findings were
observed for the CGGA, [86]GSE108474, and [87]GSE43378 (Fig. [88]2f–h,
Fig. S2a–d).
Univariate and multivariate Cox analyses of risk scores
In the TCGA dataset, univariate analysis showed the risk score was
significantly associated with OS with a hazard ratio (HR) = 4.075, 95%
confidence interval (CI) 3.108–5.344, and P < 0.001 (Fig. [89]3a). It
was also significantly associated with age (Fig. [90]3c), WHO tumor
grade (Fig. [91]3e), isocitrate dehydrogenase (IDH) mutant status
(Fig. [92]3g), and 1p/19q codeletion status (Fig. [93]3h) with each
having a P-value < 0.001. Further multivariate Cox analysis showed the
risk score was an independent factor of OS (HR = 3.330, 95% CI
1.773–6.254, P < 0.001; Fig. [94]3a). Similar to that of the TCGA
dataset, risk score was also found to be associated with OS and an
independent prognostic factor of the CGGA, [95]GSE108474, and
[96]GSE43378 datasets based on both univariate and multivariate Cox
analyses (Fig. [97]3b, Fig. S3a, b).
Fig. 3.
[98]Fig. 3
[99]Open in a new tab
Univariate and multivariate Cox regression analysis of lower-grade
glioma (LGG) prognosis in The Cancer Genome Atlas (TCGA) (a) and
Chinese Glioma Genome Atlas (CGGA) (b) datasets. Evaluation of
clinicopathological characteristics in the TCGA dataset, including age
(c), gender (d), World Health Organization (WHO) grade (e),
O-6-methylguanine-DNA methyltransferase (MGMT) methylation status (f),
isocitrate dehydrogenase (IDH) mutation status (g), and 1p/19q
codeletion status (h). *p < 0.05, **p < 0.01 and ***p < 0.001
Risk score stratified by clinicopathological characteristics
We also wanted to determine whether the risk scores correlated with
various clinicopathological characteristics of patients with LGG.
Results showed that patient age > 40, WHO grade III, unmethylated
O-6-methylguanine-DNA methyltransferase (MGMT), wild-type IDH, and
non-1p/19q codeletion were found to have higher risk scores
(Fig. [100]3c, e–h; all P < 0.001), but no significant difference was
observed in terms of risk score between females and males (P = 0.36;
Fig. [101]3d). Similar results were observed for the CGGA dataset (Fig.
S4).
Risk score benefit to current molecular classification
According to the existing molecular classification systems, patients
exhibited apparently different OS when stratified by WHO grading, IDH
mutation status, and 1p/19q codeletion status (P < 0.001, Fig. [102]4a,
c, d), but no significant difference was observed between subgroups
stratified by MGMT methylation status (P = 0.091, Fig. [103]4b).
Further analyses showed that patients in each subtype could also be
divided into two different groups in terms of OS according to their
risk score (all P < 0.05, Fig. S5a–d, f–i). Of note, patients with WHO
grade II tumors and high-risk scores had a parallel prognosis to that
of patients with WHO grade III tumors and low-risk scores
(Fig. [104]4e). A similar finding was observed in the existing
classification system of 1p/19 codeletion (Fig. [105]4h).
Fig. 4.
[106]Fig. 4
[107]Open in a new tab
Kaplan–Meier curves of multiple subgroups, including World Health
Organization (WHO) grading (a), O-6-methylguanine-DNA methyltransferase
(MGMT) methylation status (b), isocitrate dehydrogenase (IDH) mutation
status (c), and 1p/19q codeletion status (d). The risk score was
further explored to supplement WHO grading (e) and molecular subtypes,
including MGMT methylation status (f), IDH mutation status (g), and
1p/19q codeletion status (h). Association between tumor mutation burden
(TMB) and ferroptosis-related lncRNA (FRlncRNA) low-risk and high-risk
scores (i). Kaplan–Meier curves of TMB (j) and TMB with low-risk and
high-risk score (k)
TMB analysis
The TMB score in the subgroup with a high-risk score was higher than
that in the subgroup with a low-risk score (P < 0.001; Fig. [108]4i,
Table S4) and patients with high TMB were found to have worse OS than
those with low TMB (P < 0.001, Fig. [109]4j). K–M curves showed that
patients with high-risk score had a worse OS than those with low-risk
score both in the subgroup of high TMB and low TMB (both P < 0.05, Fig.
S5e, j). However, further analysis showed that patients with high TMB
and low-risk scores had a similar OS compared with that of patients
with low TMB and high-risk scores (Fig. [110]4k).
Functional analysis
A total of 739 DEGs were identified between the low-risk and high-risk
groups of the TCGA dataset and 63 DEGs between the low-risk and
high-risk groups of the CGGA dataset. As shown in Fig. [111]5a, GO
enrichment revealed many enriched terms clearly related to immunity,
including humoral immune response (biological process), MHC protein
complex (cellular component), and antigen binding (molecular function).
KEGG analysis identified immune-related pathways, such as complement
and coagulation cascades, antigen processing and presentation, and
intestinal immune network for IgA production, and DEGs accompanied by
an activated PI3K-Akt signaling pathway were also enriched
(Fig. [112]5b). The results were confirmed with similar findings for
the CGGA dataset (Fig. S6).
Fig. 5.
[113]Fig. 5
[114]Open in a new tab
The most significant gene ontology (GO) enrichment (a) and Kyoto
Encyclopedia of Genes and Genomes (KEGG) pathways (b) in The Cancer
Genome Atlas (TCGA) dataset. Multiple gene set enrichment analysis
(GSEA) for metabolism-related (c) and immune-related (d) pathways of
the same dataset
Functional annotation was performed using GSEA. DEGs were only enriched
in the high-risk group (P-value < 0.05, Table S5). Metabolism related
signaling pathways including galactose, nicotinate, phenylalanine,
starch, and sucrose were enriched in the high-risk score group
(Fig. [115]5c), as well as immunity related signaling pathways such as
antigen processing and presentation, B-cell receptor signaling pathway,
and T-cell receptor signaling pathway, among others (Fig. [116]5d).
Subsequent correlation analyses revealed that the PI3K-Akt signaling
pathway exhibited a significant positive correlation with the risk
score, ferroptosis, as well as the previously mentioned immune pathways
(e.g. B-cell receptor signaling pathway, T-cell receptor signaling
pathway, and), immune cell infiltration (e.g. macrophage cell), and
immune checkpoint expression (e.g. PD-1, PD-L1, and CTLA4) within both
the dataset of TCGA and CGGA (P < 0.05, Fig. S7, Fig. S8). Similar
results were confirmed in the external validation dataset
[117]GSE108474 (Fig. S9).
Immune response and risk score
As shown in Fig. [118]6a, many tumor-infiltrating immune cells were
found to be associated with evaluated FRlncRNA risk scores using
different infiltration estimating tools, including XCELL, TIMER,
QUANTISEQ, MCPcounter, EPIC, and CIBERSORT. In particular, estimations
of all macrophage subtypes in the high-risk group were higher than
those in the low-risk group (P < 0.001, Fig. [119]6b). Estimations of
most CD4^+ T cells in the high-risk group were also higher than those
in the low-risk group (P < 0.001, Fig. [120]6c). However, all
estimations of NK cells in the high-risk group were lower than those in
the low-risk group (P < 0.001, Fig. [121]6e). Notably, no definite
difference was observed in terms of CD8^+ T cells between the high-risk
and low-risk groups (Fig. [122]6d).
Fig. 6.
[123]Fig. 6
[124]Open in a new tab
Correlation coefficient between the tumor-infiltrating immune cells and
ferroptosis-related lncRNA (FRlncRNA) risk score (a). The boxplots show
the differences in infiltrating immune cell subtypes, including
macrophages (b), CD4+ T cells (c), CD8+ T cells (d), and natural killer
(NK) cells (e) between the low-risk score (yellow) and high-risk score
(blue) groups. Association of microenvironment score (f) and PD-1,
PD-L1, and CTLA4 expression (g) between the low-risk score (yellow) and
high-risk score (blue) groups. *p < 0.05, **p < 0.01 and ***p < 0.001
The tumor environment and immune checkpoints between the high-risk and
low-risk groups were also explored using expression data (ESTIMATE).
Results for the tumor environment showed that stromal score, immune
scores, and total score were each higher in the high-risk score group
compared to those in the low-risk score group (P < 0.001,
Fig. [125]6f). Immune checkpoints results indicated that PD-1, PD-L1,
and CTLA4 were all overexpressed in the high-risk score group compared
with those in the low-risk score group (P < 0.001, Fig. [126]6g).
Relationship between immune subtypes and risk score
The degree of immune infiltration in each sample was evaluated using
CIBERSORT algorithm, which served as the foundation for implementing
unsupervised clustering to identify two immune subtypes exhibiting
significantly distinct immune characteristics (Fig. S10a). Further
analysis revealed that cluster 1, despite exhibiting a higher risk
score and poorer prognosis (P < 0.05, Fig. S10b, c), demonstrated
elevated levels of immune checkpoint expression (P < 0.05, Fig.
S10d–f). This finding indicates that cluster 1 exhibits a stronger
correlation with the high risk score group, which comprises patients
who may demonstrate increased sensitivity to immunotherapy.
Nomogram construction based on FRlncRNAs risk scores
A nomogram was developed based on FRlncRNA risk score to further
predict OS of LGG patients (Fig. [127]7a). The concordance-index was as
high as 0.848 (0.803–0.892). Good calibrations were observed in terms
of the predicted vs. observed survival rates at 1, 3, and 5 year
(Fig. [128]7b–d). Time-dependent ROC curves showed the current nomogram
had a higher AUC for 1-, 3-, and 5-year OS than that of the other
predictors (AUC at 1-year: 0.883, Fig. [129]7e; AUC at 3-year: 0.899,
Fig. [130]7f; and AUC at 5-year: 0.826, Fig. [131]7g).
Fig. 7.
[132]Fig. 7
[133]Open in a new tab
a Nomogram for predicting survival of lower-grade glioma (LGG)
patients. b–d Calibration plots of the nomogram for predicting the
probability of overall survival (OS) at 1, 3, and 5 years. e–g Receiver
operating characteristic (ROC) curves for the ferroptosis-related
lncRNA (FRlncRNA) risk score, patient age, patient gender, tumor grade,
O-6-methylguanine-DNA methyltransferase (MGMT) methylation status,
isocitrate dehydrogenase (IDH) mutation status, and 1p/19q codeletion
status and a nomogram for predicting 1-, 3-, and 5-year OS
Discussion
Ferroptosis and lncRNAs have both been well-concerned topics over the
past decade and LGG is a highly heterogeneous molecular malignant
tumor. In the present study, we initially identified a novel prognostic
signature comprising nine FRlncRNAs to predict the long-term outcomes
of patients with LGG. Further analysis showed that the risk score was
also correlated with clinicopathological characteristics of LGG
patients. Comprehensive analyses revealed that the risk score was
highly associated with anti-tumor immune. These findings revealed that
the novel risk score could be taken as a promising prognostic biomarker
of LGG patients and ferroptosi might be a potential therapeutic target.
Overall, 1637 (10.3%) FRlncRNAs in the TCGA dataset and 493 (11.5%) in
the CGGA dataset were found to be associated with OS, indicating
FRlncRNAs might play an extensive and important role in the outcome of
LGG. The risk score exhibited better discrimination than the other
models (age, gender, grade, and status of MGMT Methtlation, IDH
mutation, 1p/19q codeletion), allowed the patients to be divided into
two subgroups that exhibited quite divergent OS both in datasets of
TCGA and CGGA (both P < 0.01). In addition, the current risk score was
also identified as an independent risk factor of OS using multivariate
Cox regression analysis both in the two datasets (both P < 0.001). All
these findings above indicated that the novel risk score could be taken
as a promising prognostic biomarker of LGG patients.
The nine FRlncRNAs included five protective ones (TMEM254-AS1,
SLC25A21-AS1, DGCR9, LINC00844, and LINC00641) and four risk ones
(AGAP2-AS1, FAM181A-AS1, CRNDE, and PAXIP1-AS2). Previous studies have
found that SLC25A21-AS1 can induce multidrug resistance in
nasopharyngeal carcinoma through miR-324-3p/IL-6 [[134]26]. Meanwhile,
DGCR9 expression is upregulated in gastric cancer and is involved in
cell proliferation and migration and affects glucose metabolism
[[135]27]. However, in the current study, these two lncRNAs were
determined to be protective factors and deserve further validation.
LINC00844 has been found to inhibit prostate cancer and hepatocellular
carcinoma proliferation by regulating the expression of GSTP1 and NDRG1
[[136]28, [137]29], respectively. Yang et al. found that LINC00641 can
suppress glioma progression through the miR-4262/NRGN axis, which
provides a novel therapeutic target [[138]30]. Meanwhile, Yan et al and
Luo et al reported that AGAP2-AS1 can promote the proliferation,
migration, and invasion of glioma cells by regulating PTN, EZH2, and
LSD1 [[139]31]. FAM181A-AS1 and CRNDE have been found to promote glioma
development by enhancing ZRANB2 expression and inhibiting Bcl-2 and
Wnt2 expression, respectively [[140]32, [141]33]. Considering there are
no relevant studies on PAXIP1-AS2 and its similar origin lncRNAs,
lncRNA PAXIP1-AS1 was found to promote glioma cell invasion and
angiogenesis by upregulating the ETS1/KIF14 axis, the prognostic value
of PAXIP1-AS2 was reliable [[142]34]. In summary, the nine FRlncRNAs we
identified may also be considered as potential novel targets for the
management of LGG.
It is well known that LGG is highly heterogeneous and the existing
molecular subtyping is inadequate to guide clinical practice [[143]2].
In the current study, we found that the risk score was quite correlated
with clinicopathological characteristics of LGG patients, including
age, WHO grading, TMB, and status of MGMT methylation status, IDH
mutation status, and 1p/19q codeletion (all P < 0.01). And then, we
found that the current risk score could also re-categorize subgroup
patients stratified by the current molecular/gene classification
system. Furthermore, patients with WHO grade II tumors and high-risk
scores were found to have a parallel prognosis with patients with WHO
grade III tumors and low-risk score; similar findings were observed in
the classification system of 1p/19q codeletion status. These findings
might explain why LGG patients of the same subtype often have different
outcomes, and the current risk score could be taken as a biomarker to
supplement the existing molecular classification system.
Ferroptosis is not only associated with tumorigenesis and aggressive
progression, but also associated with treatment response including
immune therapy [[144]4, [145]35–[146]37]. To the best of our knowledge,
all trials of immune checkpoint inhibitors to treat glioma have failed,
but the mechanism remains unclear. Similar to FRGs in previous reports
[[147]8, [148]10], the current risk score was also associated with
immunity. Immune cells including macrophages, CD8+ T cells, and NK
cells were found to be enriched in subgroup with the high-risk score
using seven different evaluation tools (all P < 0.05), which indicated
that the conclusion was robust. On the other hand, immune-related
pathways such as antigen processing and presentation, the B-cell
receptor signaling pathway, and the T-cell receptor signaling pathway
were also activated in patients with high-risk scores, which suggested
that inhibitors of these pathways may be considered as potential
therapeutic targets [[149]38]. In addition, PD-1, PD-L1, and CTLA4 were
all found in this study to be over-expressed in patients with high-risk
scores, which provides insight into immune therapy of LGG and the
screening the potential patients that may benefit from a particular
treatment, which are key for a successful trial.
This study had several limitations. First, the definition of LGG
remains controversial, even though it has been widely adopted in
bioinformatics-related studies [[150]39–[151]41]. Second, although the
data from the TCGA and CGGA databases were adjusted before subsequent
analyses, they were obtained from different batches. Third, the current
risk score exhibited good calibration for the TCGA and CGGA datasets,
but it lacked validation using external real-world data. Fourth, the
mechanism of tumor aggression is highly complex and a single hallmark
to predict OS is unlikely. Fifth, additional interesting findings
between the current risk score and the immune system were observed, but
future experiments are warranted.
Conclusion
In summary, we determined a risk score based on nine FRlncRNAs, which
could not only be taken as a promising biomarker of prognosis of LGG,
but also supplement the existing molecular classification systems. Our
findings also provide clues for future research regarding therapeutic
targets, including immune checkpoint inhibitors. However, the
underlying mechanisms require further validation using in vitro and in
vivo experiments.
Supplementary Information
[152]12672_2024_1587_MOESM1_ESM.tif^ (1.8MB, tif)
Supplementary Material 1. Fig. S1 (a) LASSO coefficient profiles of the
expression of the candidate FRlncRNAs. (b) Selection of the penalty
parameter (Lambda) in the LASSO model.
[153]12672_2024_1587_MOESM2_ESM.tif^ (4.5MB, tif)
Supplementary Material 2. Fig. S2 Kaplan–Meier curves (a), risk plot
distribution and survival status (b) of GSE108474 dataset. Kaplan–Meier
curves (c), risk plot distribution and survival status (d) of GSE43378
dataset.
[154]12672_2024_1587_MOESM3_ESM.tif^ (3.4MB, tif)
Supplementary Material 3. Fig. S3 Univariate and multivariate Cox
regression analysis of lower-grade glioma (LGG) prognosis in GSE108474
dataset (a) and GSE43378 (b) datasets.
[155]12672_2024_1587_MOESM4_ESM.tif^ (1.9MB, tif)
Supplementary Material 4. Fig. S4 Clinicopathological characteristics
evaluation including age (a), gender (b), WHO grade (c), MGMT
methylation status (d), IDH mutation status (e) and 1p/19q codeletion
status (f) in CGGA dataset. *p < 0.05, **p < 0.01 and ***p < 0.001.
[156]12672_2024_1587_MOESM5_ESM.tif^ (1.5MB, tif)
Supplementary Material 5. Fig. S5 Patients in each subtype could also
be divided into two different groups in terms of OS according to their
risk score. WHO grading (a, f), MGMT methylation status (b, g), IDH
mutation status (c, h), 1p/19q codeletion status (d, i), and TMB (e,
j).
[157]12672_2024_1587_MOESM6_ESM.tif^ (2MB, tif)
Supplementary Material 6. Fig. S6 GO enrichment (a) and KEGG pathways
(b) in the CGGA set.
[158]12672_2024_1587_MOESM7_ESM.tif^ (9.4MB, tif)
Supplementary Material 7. Fig. S7 Scatter plot and heat map showed the
correlations among PI3K-Akt signaling pathway, risk score, ferroptosis,
as well as the previously mentioned immune pathways (e.g. B-cell
receptor signaling pathway, T-cell receptor signaling pathway), immune
cell infiltration (e.g. macrophage cell and T cell CD4+), and immune
checkpoint expression (e.g. PD-1, PD-L1, and CTLA4) in TCGA cohort.
[159]12672_2024_1587_MOESM8_ESM.tif^ (9.9MB, tif)
Supplementary Material 8. Fig. S8 Scatter plot and heat map showed the
correlations among PI3K-Akt signaling pathway, risk score, ferroptosis,
as well as the previously mentioned immune pathways (e.g. B-cell
receptor signaling pathway, T-cell receptor signaling pathway), immune
cell infiltration (e.g. macrophage cell and T cell CD4+), and immune
checkpoint expression (e.g. PD-1, PD-L1, and CTLA4) in CGGA cohort.
[160]12672_2024_1587_MOESM9_ESM.tif^ (8.5MB, tif)
Supplementary Material 9. Fig. S9 Scatter plot and heat map showed the
correlations among PI3K-Akt signaling pathway, risk score, ferroptosis,
as well as the previously mentioned immune pathways (e.g. B-cell
receptor signaling pathway, T-cell receptor signaling pathway), immune
cell infiltration (e.g. macrophage cell and T cell CD4+), and immune
checkpoint expression (e.g. PD-1, PD-L1, and CTLA4) in GSE108474
cohort.
[161]12672_2024_1587_MOESM10_ESM.tif^ (5.8MB, tif)
Supplementary Material 10. Fig. S10 (a) Heatmap representation of
unsupervised clustering. (b) Kaplan-Meier curves of OS with two
clusters. (c)Boxplot shows the distribution of risk score between C1
and C2. (d)Boxplot shows the distribution of PD-1 expression between C1
and C2. (e)Boxplot shows the distribution of PD-L1 expression between
C1 and C2. (f)Boxplot shows the distribution of CTLA4 expression
between C1 and C2.
[162]12672_2024_1587_MOESM11_ESM.xlsx^ (38.1KB, xlsx)
Supplementary Material 11. Table S1 A total of 382 ferroptosis-related
genes.
[163]12672_2024_1587_MOESM12_ESM.xlsx^ (2.9MB, xlsx)
Supplementary Material 12. Table S2 A total of 1637 ferroptosis-related
lncRNAs were identified by co-expression analysis.
[164]12672_2024_1587_MOESM13_ESM.xlsx^ (28.2KB, xlsx)
Supplementary Material 13. Table S3 The 81 intersected
ferroptosis-related prognostic lncRNAs.
[165]12672_2024_1587_MOESM14_ESM.xlsx^ (31KB, xlsx)
Supplementary Material 14. Table S4 The risk score and TMB information
in TCGA cohort.
[166]12672_2024_1587_MOESM15_ESM.xlsx^ (15.7KB, xlsx)
Supplementary Material 15. Table S5 The enriched pathways of the
high-risk score groups in TCGA cohort.
[167]12672_2024_1587_MOESM16_ESM.xlsx^ (567.7KB, xlsx)
Supplementary Material 16. Table S6 The infiltration estimation for
TCGA cohort in LGG.
Acknowledgements