Advances in Clinical and Experimental Medicine

Title abbreviation: Adv Clin Exp Med
JCR Impact Factor (IF) – 2.1 (5-Year IF – 2.0)
Journal Citation Indicator (JCI) (2023) – 0.4
Scopus CiteScore – 3.7 (CiteScore Tracker 3.8)
Index Copernicus  – 171.00; MNiSW – 70 pts

ISSN 1899–5276 (print)
ISSN 2451-2680 (online)
Periodicity – monthly

Download original text (EN)

Advances in Clinical and Experimental Medicine

2024, vol. 33, nr 8, August, p. 843–856

doi: 10.17219/acem/172576

Publication type: original article

Language: English

License: Creative Commons Attribution 3.0 Unported (CC BY 3.0)

Download citation:

  • BIBTEX (JabRef, Mendeley)
  • RIS (Papers, Reference Manager, RefWorks, Zotero)

Cite as:


Zheng G, Jiang Q, Jiang C, Zhu X, Wang Y. Identification and verification on prognostic index of glioblastoma immune-related lncRNAs. Adv Clin Exp Med. 2024;33(8):843–856. doi:10.17219/acem/172576

Identification and verification on prognostic index of glioblastoma immune-related lncRNAs

Guofu Zheng1,A,B,C,D,E,F, Qingsong Jiang1,C,E, Cai Jiang1,B,C, Xingxin Zhu2,A,B,D,E, Yongjie Wang3,A,E,F

1 Department of Neurosurgery, Kecheng People’s Hospital, Quzhou, China

2 Department of Thoracic Surgery, The First Affiliated Hospital of Zhejiang University, School of Medicine, Hangzhou, China

3 Department of Neurosurgery, The Second Affiliated Hospital of Zhejiang University, School of Medicine, Hangzhou, China

Abstract

Background. Glioblastoma (GBM) is the most common cause of primary brain malignancy. Recently, many immune-related long noncoding ribonucleic acids (ir-lncRNAs) are indicated to be closely related to the regulation of the immune microenvironment and immune cell infiltration of GBM.

Objectives. Through the joint analysis of multiple public databases, key ir-lncRNAs in GBM were screened. The ir-lncRNAs were used to construct risk-scoring models and promote the development of novel GBM biomarkers.

Materials and methods. In this study, we performed a three-way Venn analysis combined with a least absolute shrinkage and selection operator (LASSO) regression analysis on all lncRNAs in The Cancer Genome Atlas (TCGA), the Chinese Glioma Genome Atlas (CGGA) and Imm-Lnc datasets, and identified 10 ir-lncRNAs. Multivariate Cox analysis was used to calculate the coefficient and construct a risk-scoring model.

Results. By plotting calibration curves and receiver operating characteristic (ROC) curves, the model showed excellent prediction results. Based on the Tumor Immune Estimation Resource (TIMER) database, the correlation analysis showed that 10 ir-lncRNAs risk scores were related to immune cell infiltration. The enrichment analysis was subsequently performed, which showed that these ir-lncRNAs played an important role in the progression of GBM. Among the 10 lncRNAs, we found that AL354993.1 was highly expressed in GBM, had not been reported, and was shown to be closely related to GBM progression.

Conclusions. In conclusion, the 10 ir-lncRNAs have the potential to predict the prognosis of GBM patients and may play a vital role in the progression of the disease.

Key words: immunity, lncRNA, prognostic model, glioblastoma, gene function enrichment analysis

Background

The most aggressive primary brain malignancy originates from oligodendrocyte or astrocyte precursor cells and is known as glioblastoma (GBM). Although accurate surgical resection, radiation and adjuvant chemotherapy are now the conventional treatments for GBM, the prognosis is still poor, and the median survival is just 8–15 months.1 Immunomodulatory therapy is a new and effective treatment option.2 The stemness features of GBM are strictly connected to immune infiltration,3 meaning neoadjuvant anti-programmed cell death protein 1 (PD-1) checkpoint blocking immunotherapy might improve the prognosis of properly selected GBM patients.4 However, GBM often exhibits severe local immunosuppression, which limits the efficacy of immunotherapy strategies.5 To further explain the mechanisms of immune regulation in GBM and offer a theoretical basis for GBM immunological treatment, we evaluated effective immune-related prognostic factors and constructed a prognostic model for GBM patients.

Long noncoding ribonucleic acids (lncRNAs) are a group of transcripts with a length of more than 200 nt that primarily function as regulators rather than protein-coding genes.6 The lncRNAs perform their biological functions in a variety of ways, including alternative splicing, transcription regulation, messenger RNA (mRNA) stability maintenance, chromatin modification, functional micropeptides, and interaction with proteins or small RNAs.7, 8, 9 The lncRNAs are also crucial for GBM progression. The lncRNA HNF1A-AS1 was shown to drive GBM progression through the microRNA (miR)-22-3p/alpha-enolase 1 (ENO1) axis.10 Indeed, lncRNA miR155HG has been shown to promote GBM progression by upregulating annexin A2 (ANXA2) as a competing endogenous RNA (ceRNA) of the tumor suppressor miR-185.11

Immune-related lncRNAs (ir-lncRNAs) are involved in regulating the GBM immune microenvironment and have unique prognostic value. According to reports, lncRNA AC003092.1 is connected to the immunosuppressive environment in GBM.12 Moreover, maternally expressed 3 (MEG3) levels are negatively associated with dendritic cell infiltration and positively correlated with infiltrating CD8+ T cells. The survival of GBM patients was also significantly correlated with the degree of MEG3 variation in copies.13 The heat shock protein family A member 7 (HSPA7) lncRNA was found to promote macrophage recruitment to the GBM tumor microenvironment and had a great prognostic value.14 However, few investigations have established prognostic models based on the identification of ir-lncRNAs in GBM.

Objectives

The ir-lncRNAs obtained from the Chinese Glioma Genome Atlas (CGGA; http://www.cgga.org.cn/) and The Cancer Genome Atlas (TCGA; https://www.cancer.gov/ccg/research/genome-sequencing/tcga) were examined. The clinical prognostic model of GBM was developed after the least absolute shrinkage and selection operator (LASSO) algorithm identified the most critical lncRNAs. Additionally, the underlying pathway of ir-lncRNAs in GBM was investigated.

Materials and methods

Data and resources

Both TCGA-GBM (n = 166) and CGGA cohorts (n = 140) were used as public transcriptome datasets in our analysis.15 Any case with a survival information null value was eliminated. The UCSC Xena database was used to retrieve the clinical information and fragments per kilobase per million (FPKM) data for the TCGA-GBM cohort (https://xenabrowser.net/). Transcripts per kilobase million (TPM) values were obtained from all FPKM data. The RNA-sequencing (RNA-seq) data of 140 specimens were retrieved from the CGGA data collection in addition to the clinical data for use as a validation set. The TCGA database provided gene mutation data (MAF files) for the TCGA-GBM group. The proportional hazards assumption test, linearity assumption test and multicollinearity test assessed the TCGA and CGCA cohorts (Supplementary Fig. 1–3).

Detection of immune-related lncRNA prognostic signature

The ImmLnc database (http://bio-bigdata.hrbmu.edu.cn/ImmReg/index.jsp) has collected 3115 GBM ir-lncRNAs.16 By evaluating the intersection of lncRNAs among the TCGA, CGGA and ImmLNC datasets, we selected ir-lncRNAs. To prevent overfitting and examine the ideal ir-lncRNA signature for estimating the overall survival of GBM individuals, the LASSO was selected. The LASSO regression analysis was performed using the “glmnet” R program. Every sample’s risk score was determined from the formula: risk score = expression value of lncRNA 1 × coefficient + expression value of lncRNA 2 × coefficient + … + expression value of lncRNA n × coefficient. Then, depending on the middle threshold of the risk score, GBM patients were allocated into elevated- and reduced-risk cohorts. The “Survival” program of R software’s area under the curve (AUC) function was employed to verify the specificity and sensitivity of the immune-related signature.

Nomogram

To anticipate the 1-, 2- and 3-year survival rates, a nomogram was developed after the independent prognostic parameters were identified. Receiver operating characteristic (ROC) curves were employed to assess the effectiveness of the model. Additionally, calibration plots were shown utilizing the rms tool to compare the model-predicted survival with the actual survival probability.

Gene set enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) and Gene Ontologygy (GO; https://geneontology.org/) mechanisms that positively related to elevated- or reduced-risk scores were investigated using gene set enrichment analysis (GSEA; https://www.gsea-msigdb.org/gsea/index.jsp). Molecular Signatures Database gene sets were obtained. Typically, 1000 permutations were used in the analysis, and pathways with a false discovery rate (FDR) of less than 0.25 were detected.

Estimation of cancer immune microenvironment

The Tumor Immune Estimation Resource (TIMER) (timer.cistrome.org/) platform17 was employed to investigate the connections between risk score and immune infiltrates, such as B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells, as demonstrated by the purity-corrected partial Spearman approach. Depending on gene expression patterns, the Estimation of STromal and Immune Cells in MAlignant Tumor tissues using Expression data (ESTIMATE) program evaluated the stromal scores and immune scores.18 In addition, we utilized the tumor immune dysfunction and exclusion (TIDE algorithm; http://tide.dfci.harvard.edu/) to evaluate each participant’s potential reaction to immune checkpoint inhibitor (ICI) treatment.19

Statistical analyses

The TCGA-GBM cohort contains 166 tumor samples from patients with GMB, while the CGGA cohort includes 140 tumor samples from patients with GMB and normal tissue from 20 patients. For the TCGA-GBM cohort, we divided patients into 2 groups based on the risk score: high-risk group (n = 83) and low-risk group (n = 83). In addition, we used the TIDE algorithm to predict the responsiveness of patients in the TCGA-GBM cohort to immunotherapy, with 48 patients evaluated as responders and 120 evaluated as non-responders. According to isocitrate dehydrogenase 1 (IDH1) mutation status, patients in the CGGA cohort were divided into 2 groups: wild-type (wt) group (n = 100) and mutant (mut) group (n = 40). Based on the co-deletion status of x1p19q, patients in the CGGA cohort were divided into 2 groups: Non_codel group (n = 128) and codel group (n = 12). Additionally, patients in the CGGA cohort were divided into 2 groups based on the expression level of AL354993.1: low-expression group (n = 70) and high-expression group (n = 70).

Continuous variables were tested for normality using Kolmogorov–Smirnov or Shapiro–Wilk tests. When the sample size was ≤50, the Shapiro–Wilk test was employed. Otherwise, the Kolmogorov–Smirnov test was used. The variables were considered to conform to a normal distribution when p > 0.05. When performing a difference analysis for 2 sets of variables, an F-test was used to evaluate the homogeneity of variance between the 2 groups. The variance between the 2 variables was considered equal when p > 0.05. The results of the tests assessing the assumptions are provided in Supplementary Tables.

When the variables conformed to a normal distribution and the variance was equal, Student’s t-test was used to compare between-group differences. If at least one of the assumptions was violated, the Mann–Whitney U test was used to compare the differences between the groups. For Student’s t-test, we also calculated the test values and degrees of freedom (df). For the Mann–Whitney U test, we calculated the U and Z values. The χ2 or Fisher’s exact tests were used to compare the differences between the groups for categorical variables. When the total sample size was >40 and the minimum theoretical frequency was >5, the χ2 test was used. When the total sample size was >40 and 5>, and the minimum theoretical frequency was >1, the corrected χ2 test was used. If the total sample size was <40 or the minimum theoretical frequency was <1, Fisher’s exact test was used. The results of tests assessing the differences between the groups are presented in tabular form and illustrated using box-and-whisker plots, which contain 5 lines representing the estimated maximum upper quartile (QU), median lower quartile (QL), and estimated minimum of the data from top to bottom. Outliers were defined as a value greater than QU + 1.5 × QR or less than QL – 1.5QR, where QR = QU – QL.

The Cox proportional hazards model was used to evaluate the impact of clinical parameters on patient survival time. In terms of parameter selection, we evaluated the clinical parameters common to both TCGA and CGGA cohorts, and age and gender were included in both. Therefore, we included 3 parameters: patient age, gender and the risk score calculated using the Cox regression model. First, we performed the proportional hazards assumption test on all 3 parameters (Supplementary Fig. 1). When the Schoenfeld individual test p-value was less than 0.05, the proportional hazards assumption was considered valid. For the continuous parameters, age and risk score, we performed a linearity assumption test (Supplementary Fig. 2) and a multicollinearity test (Supplementary Fig. 3). When the fitted curve was approximately linear, the linearity assumption was considered valid. When the p-value of the correlation was less than 0.05, these 2 parameters were considered to have no multicollinearity. For each Cox regression result, we calculated Harrell’s compliance index as the goodness-of-fit.

The sample size is labeled in the figure legends. The Kaplan–Meier technique (R package survival) was employed to create overall survival (OS) curves, and the log-rank test was used to assess alterations between the curves. This study used p < 0.05 as a statistically significant criterion. Statistical analyses employed R 3.6.2 (R Foundation for Statistical Computing, Vienna, Austria) and IBM SPSS v. 26.0 for Windows (IBM Corp., Armonk, USA).

Results

Determination of ten survival-related
ir-lncRNAs in GBM

There were 14,038 lncRNAs in the TCGA-GBM database, 1,082 lncRNAs in the CGGA database and 3,175 ir-lncRNAs in the ImmLnc database. After three-way Venn analysis, 276 ir-lncRNAs that coexisted in all 3 databases were identified (Figure 1A). To further screen for effective prognostic markers, we used LASSO regression to finally select 10 ir-lncRNAs associated with GBM prognosis, namely CHRM3-AS2, LEF1-AS1, AL354993.1, AC034243.1, LINC00452, NDUFB2-AS1, LINC00571, UNC5B-AS1, AC009093.1, and HOXC-AS2 (Figure 1B,C). The Kaplan–Meier examination indicated that elevated expression of each ir-lncRNA was substantially related to poorer OS in the TCGA database, except for NDUFB2-AS1 and LINC00571 (Figure 1C).

Progression and validation
of the ten survival-related ir-lncRNA signatures for survival anticipation

Based on the multivariate Cox regression model, the above 10 ir-LncRNAs were integrated to create a risk score model in the TCGA database. In the TCGA, the Kaplan–Meier examination revealed that OS was considerably poorer in raised-risk participants than in decreased-risk individuals (p < 0.0001, Figure 2A). The risk scores and survival status of every GBM specimen were represented by the risk curve and scatterplot, respectively. The samples in the elevated-risk group had greater risk ratings and mortality rates than those in the reduced-risk cohort (Figure 2C).

The ROC analysis revealed that the AUC for 1-, 2- and 3-year survival were 0.716, 0.803 and 0.852, respectively (Figure 2E). Therefore, the risk score model constructed using the 10 ir-lncRNAs was effective in predicting GBM prognosis. Additionally, we developed and assessed the risk score model using the CGGA database as the validation set and acquired comparable outcomes (Figure 2B,D,F).

Establishment and evaluation of nomograms by the risk scores and the prognostic value of medical variables

A univariate Cox regression examination was conducted to determine whether the risk score model for the 10 ir-lncRNAs and GBM-related clinical variables were prognostic factors. The results indicated that risk scores (p < 0.001; hazard ratio (HR) = 13.68) and patient age (p < 0.001; HR = 1.02) were strongly correlated with OS in the TCGA database (Figure 3A). The CGGA database was also used as a validation set for univariate Cox regression analysis, and the ir-lncRNA model was found to be a substantial risk factor for GBM individuals (Figure 3B).

The nomogram with age, gender and risk score was developed for the prediction of patient prognosis in the TCGA dataset. We collected survival information from all patients to anticipate the 1-, 2- and 3-year OS (Figure 3C). The calibration curve of the 1-, 2- and 3-year OS showed that the nomogram had excellent prognostic value (Figure 3E). The AUC of the 1-, 2- and 3-year OS were 0.79, 0.82 and 0.91, respectively, according to the ROC curve analysis (Figure 3F). Similar prediction outcomes were achieved in the CGGA validation set, supporting the nomogram’s effective prediction ability over the risk score model (Figure 3D,G,H).

Enrichment analysis identified
ir-lncRNA-related biological functions and signaling pathways in GBM

Based on the risk score, GBM participants in the TCGA were separated into elevated- and reduced-risk cohorts. Gene ontology and KEGG enrichment analysis were then performed on the genes that were differently expressed between the 2 groups. The outcomes of KEGG enrichment investigations revealed that the differentially expressed genes were primarily related to the activation of innate immune response, cellular responses to chemokines, collagen binding, cytokine secretion, humoral immune responses, the Janus kinase (JAK)-signal transducer and activator transcription (STAT) cascade, leukocyte migration, macrophage activation, modulation of lymphocyte activity, and T cell activation (Figure 4A). According to the outcomes of the GO enrichment analysis, these genes primarily participated in the signaling pathways for cell adhesion molecules, chemokine signaling, cytokine–cytokine receptor interactions, JAK-STAT signaling, natural killer cell-mediated cytotoxicity, nuclear factor-kappa B (NF-ĸB) signaling, nod-like receptor signaling, T helper (Th)17 cell differentiation, tumor necrosis factor (TNF) signaling, and toll-like receptor signaling (Figure 4B). The knowledge of the underlying pathways involving these ir-lncRNAs in the development and progression of GBM was deepened by these enriched biological processes and signaling mechanisms.

The ir-lncRNAs-related risk score was correlated with the GBM immune microenvironment

To additionally explore the function of the 10 ir-lncRNAs in the immune microenvironment, the TIMER database was used to examine the connection between the risk score and various infiltrating immune cell subpopulations in GBM. It was shown that the risk score was closely related to the infiltration of CD8+ T cells, neutrophils, macrophages, and myeloid dendritic cells, with correlation coefficients of 0.24, 0.25, 0.27, and 0.24, respectively (Figure 5A). It is well-known that chemokines efficiently control immune cell infiltration in malignancies. The expression of 16 immune-related chemokines recognized to be linked with GBM development was compared in elevated- and low-risk tissues to further examine the relationship between risk scores and chemokine secretion. The outcomes demonstrated that 14 chemokines were differentially expressed (chemokine ligand (CCL)2, CCL5, CCL17, CCL20, CCL22, CCR2, CCR4, CCR5, CCR6, CCR7, CXC motif chemokine ligand (CXCL)12, CXCL16, CXC motif chemokine receptor (CXCR)4, and CXCR6) (Figure 5B, Table 1 and Supplementary Table 1).

Further calculations of the GBM patients’ immunological, stromal and ESTIMATE scores in the TCGA dataset indicated that the immune score of the reduced-risk cohort was much less than that of the elevated-risk cohort (Figure 5C–E, Table 2 and Supplementary Table 2). A strong association between the risk score and commonly employed immunological checkpoints was determined using Pearson’s correlation analysis, indicating that the risk score may have implications for immunotherapy (Figure 5F). To further anticipate the reaction of GBM to immunotherapy, we used the TIDE algorithm. A high TIDE score means that the likelihood of responding to immunotherapy is low. According to the findings, individuals with elevated-risk GBM had reduced TIDE scores and were more likely to react to immunotherapy (Figure 5G,H, Table 2, Table 3 and Supplementary Table 2). The previous findings demonstrated that the risk score associated with ir-lncRNAs was strongly connected to the GBM immune microenvironment, and it was anticipated that this relationship might be helpful in targeting future GBM immunotherapy.

Alterations in the mutation landscape between elevated- and low-risk groups

Tumor somatic mutation of commonly mutated genes was profiled in GBM elevated-risk and reduced-risk cohorts. Most mutated genes, including TP53, PTEN, TTN, EGFR, MUC16, SPTA1, NF1, RYR2, and FLG, had high rates of somatic mutations in both raised- and reduced-risk groups. In particular, somatic mutation rates for RB1, HYDIN, ATRX, and IDH1 were greater in the low-risk cohort (Figure 6A,B). Notably, ATRX and IDH1 did not have somatic mutations in the raised-risk cohort in the TCGA dataset (Figure 6C,D). In the CGGA dataset, the patterns of risk scores for IDH1 mutation and 1p/19q codeletion stats were identified. According to the findings, those with IDH1 mutations had considerably lower risk scores than individuals with IDH1 wild type (Figure 6E, Table 4 and Supplementary Table 4). Participants who had a 1p/19q codeletion also had a lower risk score than patients who did not have a codeletion (Figure 6F, Table 4 and Supplementary Table 5).

AL354993.1 was recognized as a potential marker for GBM

The expression of 10 ir-lncRNAs was compared between normal and tumor samples in the CGGA dataset. The LEF1-AS1, HOXC-AS2 and AL354993.1 were found to be significantly overexpressed in GBM (p < 0.0001), among which AL354993.1 has not been reported (Figure 7A, Table 5 and Supplementary Table 6). The overexpression of AL354993.1 suggested a worse GBM prognosis (p = 0.029) (Figure 7B). Furthermore, the expression of AL354993.1 was lower in GBM tissues with IDH1 but not with 1p19q mutations (Figure 7C,D, Table 6 and Supplementary Table 7,8). To discover the biological function of AL354993.1, GBM subjects in the CGGA dataset were allocated into 2 cohorts with raised and decreased AL354993.1 expression, and a differential gene analysis combined with GSEA enrichment was conducted. The outcomes indicated that raised AL354993.1 expression might indicate the activation of the JAK-STAT pathway, the extracellular matrix (ECM) receptor interaction pathway, chemokine signaling mechanisms, and cytokine–cytokine receptor interaction mechanisms (Figure 7E). Hence, the relationship between AL354993.1 expression and chemokine secretion were further explored. The expression of 16 chemokines known to be associated with GBM progression was examined in reduced- and elevated-expression groups. The results showed that 13 chemokines were differentially expressed (CCL2, CCL5, CCL17, CCL20, CCR2, CCR4, CCR5, CCR6, CCR7, CXCL12, CXCL16, CXCR4, CXCR6) (Figure 7F, Table 7 and Supplementary Table 9). These outcomes indicate that AL354993.1 is related to GBM progression and might be a potential GBM biomarker.

Discussion

Glioblastoma multiforme is one of the most common primary brain tumors, and survival rates are low despite extensive treatment options (surgical resection, radiotherapy and adjuvant chemotherapy).20 In recent years, many studies carried out on immunotherapy for GBM have failed to achieve ideal results due to the heterogeneity of GBM and the cancer immunosuppressive microenvironment.21, 22 Numerous investigations have verified that lncRNAs have a function in the modulation of GBM. Several ir-lncRNAs have also been shown to be strongly linked to the modulation of the immunological milieu in GBM and the infiltration of immune cells.7, 16 Therefore, we created a medical prediction model using the TCGA and CGGA datasets and evaluated the lncRNAs most relevant to GBM to further examine the function of lncRNAs in GBM prognosis.

In this investigation, we conducted a three-way Venn examination of all lncRNAs in the TCGA, CGGA and ImmLnc datasets and found that 276 lncRNAs coexisted. The LASSO regression analysis further identified 10 key ir-lncRNAs, with the Kaplan–Meier analysis suggesting that 8 of the 10 were prognostic risk factors for GBM (CHRM3-AS2, LEF1-AS1, AL354993.1, AC034243.1, LINC00452, UNC5B-AS1, AC009093.1, and HOXC-AS2). Among these, LEF1-AS1 was found to be related to the malignant growth of GBM through increased EN2 by mir-543 via sponge absorption.23 Besides, CHRM3-AS2, an immune-associated lncRNA, has been recognized as a prognostic risk factor for ovarian cancer and cholangiocarcinoma.24, 25 The LINC00452 has been shown to promote ovarian cancer by inhibiting ubiquitin-mediated degradation through sponge absorption of mir-501-3p, thereby increasing ROCK1.26 Meanwhile, UNC5B-AS1 is involved in liver malignancy, prostate tumors and colorectal malignancy.27, 28 By mixing with HOXC13, HOXC-AS2 regulates the growth, cell death and migration of non-small cell lung tumors.29 The roles of these ir-lncRNAs in GBM are worthy of further investigation.

The coefficient was determined by employing multivariate Cox analysis, and a risk score model was developed to categorize GBM individuals into high- and reduced-risk cohorts. We discovered that the elevated-risk cohort had a lower OS rate and that the risk score model was beneficial in anticipating GBM prognosis. We also developed a nomogram that included sex, age and risk score. The nomogram demonstrated excellent predictive performance when calibration and ROC curves were plotted. Since these ir-lncRNAs have independent prognostic value in GBM OS, the risk prediction model constructed by these ir-lncRNAs also showed strong prognostic value. Therefore, we believe that monitoring these ir-lncRNAs by liquid biopsy or tumor tissue collection for GBM risk and prognosis assessment is a promising strategy.

Functional enrichment analysis was conducted to investigate the probable function and role of the 10 lncRNAs in GBM formation. They were shown to be connected with immune-related activities, including chemokine-related mechanisms, the JAK-STAT signaling pathway, cytokine-related mechanisms, and immune cell activation-related pathways, according to GO and KEGG enrichment assessment. Additionally, they were connected to the stimulation of immune-related mechanisms, such as the Toll-like receptor pathway, the TNF signaling pathway and NF-ĸB signaling. To further explore the potential impact of the 10 ir-lncRNAs on GBM immunity, we examined the relationship between risk score and immune cell infiltration in GBM tissues. The outcomes demonstrated that CD8+ T cells, neutrophils, macrophages and myeloid dendritic cells were substantially positively correlated with the risk score. The relationship between risk scores and chemokines was also investigated, with the raised-risk cohort related to a higher chemokine expression.

Immunological checkpoint blocking has become one of the frontiers of cancer immunotherapy,30 and we discovered that the risk score was substantially associated with immunological checkpoint expression. The ESTIMATE and TIDE scores also verified the relationship between risk score and immune microenvironment. Mutations in the IDH1 gene are commonly found in GBM, with the product of this mutated enzyme having a novel ability to catalyze the production of 2-hydroxyglutarate, while IDH1 mutations are shown to be associated with more favorable survival outcomes.31 In addition, 1p/19q codeletion suggested a better prognosis.32 We found higher rates of IDH1 mutation and 1p/19q codeletion in GBM individuals with reduced risk scores, while no mutations were found in the raised-risk cohort. These results may be the reason underlying the prognostic value of risk scores. Furthermore, investigation is still required to fully understand the precise pathway of the 10 ir-lncRNAs in GBM.

Among the 10 lncRNAs, we found that only 2 were significantly highly expressed in GBM compared with normal controls, namely LEF1-AS1 and AL354993.1, of which AL354993.1 has not been reported. We found that AL354993.1 had significant prognostic value in GBM and was also associated with the tumor immune microenvironment. In addition, AL354993.1 was closely related to IDH1 mutation and other clinicopathological factors. Gene enrichment results suggested that AL354993.1 may promote GBM progression through the JAK-STAT pathway, the ECM interaction pathway, chemokine pathways, and cytokine receptor binding pathways.

Limitations

There were some restrictions to this analysis since it was retrospective and dependent on open-access databases. Thus, more prospective clinical data are required for validation and future clinical applications. On the other hand, the effect of AL354993.1 was only analyzed using a database and should be verified by a series of wet experiments in the future.

Conclusions

Through the joint analysis of multiple public databases, we screened 10 key ir-lncRNAs in GBM and used them to construct risk-scoring models and promote the development of novel GBM biomarkers. In conclusion, the 10 ir-lncRNAs have the potential to predict the GBM prognosis and may play a vital role in the progression of GBM.

Supplementary data

The supplementary materials are available at https://doi.org/10.5281/zenodo.8311458. The package includes the following files:

Supplementary Fig. 1. Proportional hazards assumption test for age, gender and risk score of TCGA-GBM (A) and CGGA (B) cohort.

Supplementary Fig. 2. A,B. Linearity assumption test for age of TCGA-GBM and CGGA cohort; C,D. Linearity assumption test for the risk score of TCGA-GBM and CGGA cohort.

Supplementary Fig. 3. Multicollinearity test for age and risk score of TCGA-GBM (A) and CGGA (B) cohort.

Supplementary Table 1. The expression of chemokines between low- and high-risk groups (Figure 5B).

Supplementary Table 2. The difference of immune score, stromal score, ESTIMATE score, and TIDE score between high- and low-risk groups (Figure 5C–E,G).

Supplementary Table 3. Comparison of the risk scores between immunotherapy responders (n = 48) and non-responders (n = 118) predicted with TIDE (Figure 5H).

Supplementary Table 4. The distribution of risk score between IDH1-wild type (n = 100) and IDH1-mutation (n = 40) patients in the CGGA dataset (Figure 6E).

Supplementary Table 5. The distribution of risk score between 1p/19q-codeleted and 1p/19q-noncodeleted patients in the CGGA dataset (Figure 6F).

Supplementary Table 6. The expression of immune-related lncRNAs between normal and tumor samples in the CGGA dataset (Figure 7A).

Supplementary Table 7. The expression of AL354993.1 between IDH1-wild type and IDH1-mutation patients in the CGGA dataset (Figure 7C).

Supplementary Table 8. The expression of AL354993.1 between 1p/19q-codeleted (n = 12) and 1p/19q-noncodeleted (n = 128) patients in the CGGA dataset (Figure 7D).

Supplementary Table 9. The expression of chemokines between low and high AL354993.1 expression group (Figure 7F).

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Consent for publication

Not applicable.

Tables


Table 1. The expression of chemokines between low- and high-risk groups

Gene

Gene expression

p-value

low risk (n = 83)

high risk (n = 83)

CCL2

5.84 ±1.66

6.96 ±1.46

9.00E-06

CCL5

3.33 ±1.02

4.06 ±1.01

9.00E-06

CCL17

0.83 ±0.60

1.10 ±0.95

2.80E-02

CCL20

1.58 ±1.25

2.34 ±1.73

1.40E-03

CCL22

0.54 ±0.39

0.74 ±0.63

1.40E-02

CCR2

0.74 ±0.57

1.44 ±0.97

7.34E-08

CCR4

0.29 ±0.22

0.59 ±0.54

7.00E-06

CCR5

2.31 ±0.97

2.92 ±0.79

1.60E-05

CCR6

0.02 ±0.02

0.03 ±0.04

8.16E-04

CCR7

0.73 ±0.44

1.17 ±0.70

4.00E-06

CX3CL1

5.35 ±0.83

5.31 ±0.72

7.28E-01

CX3CR1

5.45 ±1.41

5.76 ±1.21

1.28E-01

CXCL12

3.24 ±1.10

3.82 ±1.19

1.11E-03

CXCL16

6.60 ±0.86

6.94 ±0.57

2.97E-03

CXCR4

6.48 ±0.97

6.74 ±0.75

4.79E-02

CXCR6

0.59 ±0.50

1.03 ±0.65

3.00E-06

Values in bold are statistically significant.
Table 2. The difference of immune score, stromal score, ESTIMATE score, and TIDE score

Score items

Value

p-value

low risk (n = 83)

high risk (n = 83)

Immune score

–318.36 ±453.26

113.18 ±556.25

6.00E-06

Stromal score

499.02 ±700.96

993.86 ±659.55

1.58E-07

ESTIMATE score

180.66 ±1107.5

1107.04 ±1171.6

5.00E-07

TIDE score

1.06 ±0.46

0.09 ±0.83

3.33E-16

ESTIMATE – Estimation of STromal and Immune Cells in Malignant Tumor tissues using Expression data; TIDE – tumor immune dysfunction and exclusion.
Table 3. Comparison of the risk scores between immunotherapy responders (n = 48) and non-responders (n = 118) predicted with TIDE

Score items

Status

Risk score

p-value

TIDE score

non-responder (n = 118)

0.44 ±0.23

1.10E-02

responder (n = 48)

0.61 ±0.22

TIDE – tumor immune dysfunction and exclusion.
Table 4. Different risk scores divided by IDH1 mutation
and 1p/19q codeletion

Mutation

Status

Risk score

p-value

IDH1 mutation

wt (n = 100)

0.45 ±0.24

8.29E-09

mut (n = 40)

0.17 ±0.22

1p/19q codeletion

Non_Codel (n = 128)

0.38 ±0.26

2.23E-02

Codel (n = 12)

0.20 ±0.24

wt – wild type; mut – mutant. 
Table 5. The expression of immune-related lncRNAs between normal and tumor samples in the CGGA dataset

Gene

Gene expression

p-value

normal (n = 20)

tumor (n = 140)

lEF1-AS1

0.126 ±0.131

1.008 ±1.196

3.20E-14

HOXC-AS2

0.011 ±0.036

0.335 ±0.630

1.63E-08

AL354993.1

0.601 ±0.463

1.717 ±2.179

5.86E-07

NDUFB2-AS1

0.491 ±0.171

0.570 ±0.320

9.70E-01

AC034243.1

0.372 ±0.683

0.438 ±0.726

6.99E-01

UNC5B-AS1

0.877 ±0.567

0.529 ±0.693

4.82E-02

LINC00571

0.158 ±0.253

0.238 ±0.379

2.37E-01

CHRM3-AS2

0.062 ±0.060

0.086 ±0.111

1.44E-01

LINC00452

0.020 ±0.025

0.094 ±0.247

4.17E-02

AC009093.1

0.037 ±0.030

0.052 ±0.065

7.99E-02

lncRNAs – long noncoding ribonucleic acids; CGGA – Chinese Glioma Genome Atlas. Values in bold are statistically significant.
Table 6. The different expression of AL354993.1 divided by IDH1 mutation and 1p/19q codeletion

Mutation

Status

Risk score

p-value

IDH1 mutation

wt (n = 100)

1.32 ±0.0.96

3.00E-05

mut (n = 40)

0.72 ±0.62

1p/19q codeletion

Non_Codel (n = 128)

1.17 ±0.91

1.93E-01

Codel (n = 12)

0.99 ±0.98

wt – wild type; mut – mutant.
Table 7. The expression of chemokines between low- and high- AL354993.1 expression groups

Gene

Gene expression

p-value

low expression (n = 70)

high expression (n = 70)

cCL2

5.46 ±2.11

6.82 ±1.87

9.74E-05

CCL5

2.35 ±1.10

3.18 ±1.15

2.58E-05

CCL17

0.21 ±0.27

0.28 ±0.35

1.56E-01

CCL20

0.91 ±1.29

1.70 ±1.95

5.45E-03

CCL22

0.20 ±0.19

0.24 ±0.21

2.49E-01

CCR2

0.49 ±0.47

0.98 ±0.69

2.76E-06

CCR4

0.16 ±0.19

0.40 ±0.42

1.95E-05

CCR5

1.28 ±0.78

1.89 ±0.84

1.31E-05

CCR6

0.21 ±0.22

0.34 ±0.33

4.99E-03

CCR7

0.43 ±0.36

0.80 ±0.51

2.94E-06

CX3CL1

4.44 ±0.79

4.20 ±0.71

5.55E-02

CX3CR1

3.94 ±1.49

3.49 ±1.57

8.56E-02

CXCL12

3.28 ±0.99

3.86 ±1.24

2.61E-03

CXCL16

5.13 ±1.03

5.57 ±0.96

9.31E-03

CXCR4

4.58 ±1.16

5.48 ±1.09

6.21E-06

CXCR6

0.49 ±0.48

0.95 ±0.72

2.17E-05

Values in bold are statistically significant.

Figures


Fig. 1. Identification of 10 survival-related immune long noncoding ribonucleic acids (lncRNAs) in glioblastoma. A. Venn diagram illustrates that 276 lncRNAs were shared among 3 datasets; B. Distribution of least absolute shrinkage and selection operator (LASSO) coefficients for 276 immune-related lncRNAs; C. Partial likelihood deviation of the LASSO coefficient distribution. Vertical dashed lines indicate lambda.min and lambda.lse; D. The Kaplan–Meier survival curves comparing overall (OS) survival between high- and low-expression groups of the selected immune-related lncRNAs in The Cancer Genome Atlas (TCGA) dataset
CGGA – Chinese Glioma Genome Atlas.
Fig. 2. Development and validation of 10 survival-related immune long noncoding ribonucleic acids (lncRNA) signatures for survival prediction. The Kaplan–Meier survival analysis of high-risk and low-risk patients divided by the medium value in The Cancer Genome Atlas (TCGA) (A) and Chinese Glioma Genome Atlas (CGGA) (B) datasets. The distribution of risk scores, survival time and status of patients in TCGA (C) and CGGA (D) datasets. Receiver operating characteristic (ROC) analysis of the prognostic signatures to predict 1-, 2- and 3-year overall survival (OS) in TCGA (E) and CGGA (F) datasets
Fig. 3. Establishment and evaluation of nomograms using risk scores and the prognostic value of clinical variables. A. Univariate Cox analysis showed that risk scores and age were significantly related to overall survival (OS) in The Cancer Genome Atlas (TCGA) dataset using Harrell’s compliance index (C-index) as the goodness-of-fit. For the TCGA cohort, the C-index was 0.69 (95% confidence interval (95% CI): of 0.64–0.73), and for the Chinese Glioma Genome Atlas (CGGA cohort), the C-index was 0.6 (95% CI: 0.54–0.65); B. Univariate Cox analysis showed that risk scores and age were significantly related to OS in the CGGA dataset; C. Development of the nomograms for the prediction of patient prognosis in the TCGA dataset; D. Development of the nomograms for the prediction of patient prognosis in the CGGA dataset; E. The calibration curve for 1-, 2- and 3-year OS of the nomogram in the TCGA dataset; F. Receiver operating characteristic (ROC) curves displayed the area under the curve (AUC) for 1-, 2- and 3-year OS in the TCGA dataset; G. The calibration curve for 1-, 2- and 3-year OS of the nomogram in the CGGA dataset; H. ROC curves displayed the AUC of 1-, 2- and 3-year OS in the CGGA dataset
Fig. 4. Gene set enrichment analysis between high-risk and low-risk groups in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) datasets
Fig. 5. The immune risk score was correlated with the tumor immune microenvironment. A. The correlation between risk score and immune cell infiltration; B. The expression of chemokines between low- (n = 83) and high-risk groups (n = 83) (parametric test and the Mann–Whitney U test). The difference in immune score (parametric test) (C), stromal score (parametric test) (D), and Estimation of STromal and Immune Cells in Malignant Tumor tissues using Expression data (ESTIMATE) score (Mann–Whitney U test) (E) between low- (n = 83) and high-risk groups (n = 83); F. Assessment of associations between risk scores and commonly used immunological checkpoints; G. Tumor immune dysfunction and exclusion (TIDE) scores between high- (n = 83) and low-risk group (n = 83); H. Comparison of the risk scores between immunotherapy responders (n = 48) and non-responders (n = 118) predicted with TIDE (parametric test)
Fig. 6. Differences in the mutation landscape between high- and low-risk groups. A. Oncoplots of tumor somatic mutations of frequently mutated genes in the low-risk group; B. Oncoplots of tumor somatic mutations of frequently mutated genes in the high-risk group. Mutation rates of IDH1 (C) and ATRX (D) were compared between low- and high-risk groups in The Cancer Genome Atlas (TCGA) dataset; E. The distribution of risk score between IDH1-wild type (wt) (n = 100) and IDH1-mutation (Mut) (n = 40) patients in the Chinese Glioma Genome Atlas (CGGA dataset) (Mann–Whitney U test); F. The distribution of risk score between 1p/19q-codeleted (n = 12) and 1p/19q-noncodeleted (n = 128) patients in the CGGA dataset (Mann–Whitney U test)
Fig. 7. Identification of AL354993.1 as a potential biomarker for glioblastoma. A. The expression of immune-related long noncoding ribonucleic acids (ir-lncRNAs) between normal (n = 20) and tumor samples (n = 140) in the Chinese Glioma Genome Atlas (CGGA) dataset (parametric test and the Mann–Whitney U test); B. Survival analysis of patients divided by low and high AL354993.1 expression in the CGGA dataset; C. The expression of AL354993.1 between IDH1-wild type (wt) (n = 100) and IDH1-mutation (Mut) (n = 40) patients in the CGGA dataset (Mann–Whitney U test). D. The expression of AL354993.1 between 1p/19q-codeleted (n = 12) and 1p/19q-noncodeleted (n = 128) patients in the CGGA dataset (Mann–Whitney U test); E. Gene set enrichment analysis between low and high AL354993.1 expression group in the CGGA dataset; F. The expression of chemokines between low (n = 70) and high (n = 70) AL354993.1 expression group (parametric test and Mann–Whitney U test)

References (32)

  1. Aldape K, Zadeh G, Mansouri S, Reifenberger G, Von Deimling A. Glioblastoma: Pathology, molecular mechanisms and markers. Acta Neuropathol. 2015;129(6):829–848. doi:10.1007/s00401-015-1432-1
  2. Himes BT, Geiger PA, Ayasoufi K, Bhargav AG, Brown DA, Parney IF. Immunosuppression in glioblastoma: Current understanding and therapeutic implications. Front Oncol. 2021;11:770561. doi:10.3389/fonc.2021.770561
  3. Warrier NM, Agarwal P, Kumar P. Integrative analysis to identify genes associated with stemness and imamune infiltration in glioblastoma. Cells. 2021;10(10):2765. doi:10.3390/cells10102765
  4. Sun L, Lai TJ, Prins RM. Is there a role for neoadjuvant anti-PD-1 therapies in glioma? Curr Opin Neurol. 2021;34(6):834–839. doi:10.1097/WCO.0000000000000992
  5. Ayasoufi K, Pfaller CK, Evgin L, et al. Brain cancer induces systemic immunosuppression through release of non-steroid soluble mediators. Brain. 2020;143(12):3629–3652. doi:10.1093/brain/awaa343
  6. Gandhi M, Groß M, Holler JM, et al. The lncRNA lincNMR regulates nucleotide metabolism via a YBX1–RRM2 axis in cancer. Nat Commun. 2020;11(1):3214. doi:10.1038/s41467-020-17007-9
  7. Jin X, Ge LP, Li DQ, et al. LncRNA TROJAN promotes proliferation and resistance to CDK4/6 inhibitor via CDK2 transcriptional activation in ER+ breast cancer. Mol Cancer. 2020;19(1):87. doi:10.1186/s12943-020-01210-9
  8. Liu J, Liu ZX, Wu QN, et al. Long noncoding RNA AGPG regulates PFKFB3-mediated tumor glycolytic reprogramming. Nat Commun. 2020;11(1):1507. doi:10.1038/s41467-020-15112-3
  9. Zheng S, Yang L, Zou Y, et al. Long non-coding RNA HUMT hypomethylation promotes lymphangiogenesis and metastasis via activating FOXK1 transcription in triple-negative breast cancer. J Hematol Oncol. 2020;13(1):17. doi:10.1186/s13045-020-00852-y
  10. Ma C, Wang H, Zong G, et al. EGR1 modulated LncRNA HNF1A-AS1 drives glioblastoma progression via miR-22-3p/ENO1 axis. Cell Death Discov. 2021;7(1):350. doi:10.1038/s41420-021-00734-3
  11. Wu W, Yu T, Wu Y, Tian W, Zhang J, Wang Y. The miR155HG/miR-185/ANXA2 loop contributes to glioblastoma growth and progression. J Exp Clin Cancer Res. 2019;38(1):133. doi:10.1186/s13046-019-1132-0
  12. Guo XY, Zhong S, Wang ZN, et al. Immunogenomic profiling demonstrate AC003092.1 as an immune-related eRNA in glioblastoma multiforme. Front Genet. 2021;12:633812. doi:10.3389/fgene.2021.633812
  13. Xu X, Zhong Z, Shao Y, Yi Y. Prognostic value of MEG3 and its correlation with immune infiltrates in gliomas. Front Genet. 2021;12:679097. doi:10.3389/fgene.2021.679097
  14. Zhao R, Li B, Zhang S, et al. The N6-methyladenosine-modified pseudogene HSPA7 correlates with the tumor microenvironment and predicts the response to immune checkpoint therapy in glioblastoma. Front Immunol. 2021;12:653711. doi:10.3389/fimmu.2021.653711
  15. Zhao Z, Zhang KN, Wang Q, et al. Chinese Glioma Genome Atlas (CGGA): A comprehensive resource with functional genomic data from Chinese glioma patients. Genomics Proteomics Bioinformatics. 2021;19(1):1–12. doi:10.1016/j.gpb.2020.10.005
  16. Li Y, Jiang T, Zhou W, et al. Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nat Commun. 2020;11(1):1000. doi:10.1038/s41467-020-14802-2
  17. Li T, Fu J, Zeng Z, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–W514. doi:10.1093/nar/gkaa407
  18. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4(1):2612. doi:10.1038/ncomms3612
  19. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–1558. doi:10.1038/s41591-018-0136-1
  20. Badillo-Martinez D, Nicotera N, Bodnar RJ. Onset of pain threshold changes induced by neonatal monosodium glutamate. Int J Neurosci. 1984;24(3–4):275–279. doi:10.3109/00207458409089816
  21. Patel AP, Tirosh I, Trombetta JJ, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344(6190):1396–1401. doi:10.1126/science.1254257
  22. Moserle L, Casanovas O. Anti-angiogenesis and metastasis: A tumour and stromal cell alliance. J Intern Med. 2013;273(2):128–137. doi:10.1111/joim.12018
  23. Zeng S, Zhou C, Yang DH, et al. LEF1-AS1 is implicated in the malignant development of glioblastoma via sponging miR-543 to upregulate EN2. Brain Res. 2020;1736:146781. doi:10.1016/j.brainres.2020.146781
  24. Li Y, Huo FF, Wen YY, Jiang M. Screening and identification of an immune-associated lncRNA prognostic signature in ovarian carcinoma: Evidence from bioinformatic analysis. Biomed Res Int. 2021;2021:6680036. doi:10.1155/2021/6680036
  25. Liu YJ, Hounye AH, Wang Z, Liu X, Yi J, Qi M. Identification and validation of three autophagy-related long noncoding RNAs as prognostic signature in cholangiocarcinoma. Front Oncol. 2021;11:780601. doi:10.3389/fonc.2021.780601
  26. Yang J, Wang WG, Zhang KQ. LINC00452 promotes ovarian carcinogenesis through increasing ROCK1 by sponging miR-501-3p and suppressing ubiquitin-mediated degradation. Aging. 2020;12(21):21129–21146. doi:10.18632/aging.103758
  27. Wilkie ND. The search for the practical. J Prosthet Dent. 1974;32(3):251–255. doi:10.1016/0022-3913(74)90026-2
  28. Tan SF, Ni JX, Xiong H. LncRNA UNC5B-AS1 promotes malignant progression of prostate cancer by competitive binding to caspase-9. Eur Rev Med Pharmacol Sci. 2020;24(5):2271–2280. doi:10.26355/eurrev_202003_20493
  29. Liu B, Li J, Li JM, Liu GY, Wang YS. HOXC-AS2 mediates the proliferation, apoptosis, and migration of non-small cell lung cancer by combining with HOXC13 gene. Cell Cycle. 2021;20(2):236–246. doi:10.1080/15384101.2020.1868161
  30. Galluzzi L, Humeau J, Buqué A, Zitvogel L, Kroemer G. Immunostimulation with chemotherapy in the era of immune checkpoint inhibitors. Nat Rev Clin Oncol. 2020;17(12):725–741. doi:10.1038/s41571-020-0413-z
  31. Seyfried TN, Shivane AG, Kalamian M, Maroon JC, Mukherjee P, Zuccoli G. Ketogenic metabolic therapy, without chemo or radiation, for the long-term management of IDH1-mutant glioblastoma: An 80-month follow-up case report. Front Nutr. 2021;8:682243. doi:10.3389/fnut.2021.682243
  32. Kebir S, Weber M, Lazaridis L, et al. Hybrid 11C-MET PET/MRI combined with “machine learning” in glioma diagnosis according to the revised glioma WHO classification 2016. Clin Nucl Med. 2019;44(3):214–220. doi:10.1097/RLU.0000000000002398