Wanli Yu
Yanan Ma3†
Fang Wang
Pengfei Wu
Guohua Zhang- 1Department of Neurosurgery, Gaoxin Hospital of The First Affiliated Hospital of Nanchang University, Nanchang, China
- 2Department of Neurosurgery, The First Affiliated Hospital of Nanchang University, Nanchang, China
- 3Laboratory of Medical Genetics, Harbin Medical University, Harbin, China
- 4Department of Urology, The Second Affiliated Hospital of Harbin Medical University, Harbin, China
- 5Department of Neurosurgery, The Second Affiliated Hospital of Harbin Medical University, Harbin, China
- 6The Laboratory of Artificial Intelligence and Bigdata in Ophthalmology, The Affiliated Eye Hospital of Nanjing Medical University, Nanjing, China
- 7Oncology Department, Gaoxin Hospital of The First Affiliated Hospital of Nanchang University, Nanchang, China
- 8Department of Neurosurgery, The First Affiliated Hospital of USTC, Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, China
- 9Anhui Provincial Stereotactic Neurosurgical Institute, Hefei, China
- 10Anhui Province Key Laboratory of Brain Function and Brain Disease, Hefei, China
- 11Anhui Provincial Clinical Research Center for Neurosurgical Disease, Hefei, China
- 12Central Laboratory, Gaoxin Hospital of The First Affiliated Hospital of Nanchang University, Nanchang, China
Background: Glioblastoma multiforme (GBM) is extensively genetically and transcriptionally heterogeneous, which poses challenges for classification and management. Long noncoding RNAs (lncRNAs) play a critical role in the development and progression of GBM, especially in tumor-associated immune processes. Therefore, it is necessary to develop an immune-related lncRNAs (irlncRNAs) signature.
Methods: Univariate and multivariate Cox regression analyses were utilized to construct a prognostic model. GBM-specific CeRNA and PPI network was constructed to predict lncRNAs targets and evaluate the interactions of immune mRNAs translated proteins. GO and KEGG pathway analyses were used to show the biological functions and pathways of CeRNA network-related immunity genes. Consensus Cluster Plus analysis was used for GBM gene clustering. Then, we evaluated GBM subtype-specific prognostic values, clinical characteristics, genes and pathways, immune infiltration access single cell RNA-seq data, and chemotherapeutics efficacy. The hub genes were finally validated.
Results: A total of 17 prognostically related irlncRNAs were screened to build a prognostic model signature based on six key irlncRNAs. Based on GBM-specific CeRNAs and enrichment analysis, PLAU was predicted as a target of lncRNA-H19 and mainly enriched in the malignant related pathways. GBM subtype-A displayed the most favorable prognosis, high proportion of genes (IDH1, ATRX, and EGFR) mutation, chemoradiotherapy, and low risk and was characterized by low expression of four high-risk lncRNAs (H19, HOTAIRM1, AGAP2-AS1, and AC002456.1) and one mRNA KRT8. GSs with poor survival were mainly infiltrated by mesenchymal stem cells (MSCs) and astrocyte, and were more sensitive to gefitinib and roscovitine. Among GSs, three hub genes KRT8, NGFR, and TCEA3, were screened and validated to potentially play feasible oncogenic roles in GBM.
Conclusion: Construction of lncRNAs risk model and identification of GBM subtypes based on 17 irlncRNAs, which suggesting that irlncRNAs had the promising potential for clinical immunotherapy of GBM.
Introduction
Glioblastoma multiforme (GBM) is the most frequent intracranial primary malignancy in adults. Despite standard treatment, the median survival of GBM patients is less than 14 months (1). In the latest glioma classification, molecular features are considered as classifiers in conjunction with histopathological appearance (2). Emerging biosomics studies have improved the diagnosis and treatment strategies for GBM to some extent but have not yet achieved satisfactory results due to the complex pathogenesis and molecular heterogeneity of GBM. Therefore, more studies are urgently needed to explore the mechanisms involved and to identify novel biomarkers to predict the prognosis and therapeutic effects of GBM.
Long noncoding RNA (lncRNA) is a noncoding RNA with a length of more than 200 nucleotides (3). The discovery of lncRNAs has uncovered new horizons in the pathological processes of multiple diseases, including cancer initiation and progression (4). Recent studies have shown that lncRNAs can influence the tumor immune microenvironment (TIM) by regulating inflammation and participating in immune gene expression (5, 6). For example, lncRNA nuclear-enriched abundant transcript 1 (NEAT1) affects cytokine response and induces IGs expression through the regulation of interleukin (IL)-8 transcription (7). LncRNA-Cox2 participates in inflammatory gene expression in macrophages via regulating chromatin complex remodeling (8). Zhao et al. showed that the lncRNA SNHG14/miR-5590-3p/ZEB1 positive feedback loop can regulate the PD-1/PD-L1 checkpoint to promote diffuse large B cell lymphoma progression and immune evasion (9). Increasing studies reporting on the mechanism of irlncRNAs in multiple cancers (10), the ambiguous relationship between lncRNAs, and the tumor immune microenvironment have been gradually unveiled. However, the relationship between lncRNAs and tumor immune microenvironment is rarely studied in GBM. Therefore, identification of the irlncRNAs signature may provide a new insight for predicting prognosis and individualized treatment of GBM.
In this study, we identified six key irlncRNA signatures (H19, ST3GAL6-AS1, AL162231.2, SOX21-AS1, AC006213.5, and AC002456.1), which concluded that the risk model indeed had a good predictive outcome. GBM-specific CeRNAs were constructed to predict irlncRNAs targets. GO and KEGG pathway enrichment analysis was used to explore target functions. The PPI network was performed to identify the interactions of proteins translated from mRNAs in the CeRNA network. Furthermore, GS-A showed better prognosis among the identified four GSs (A-D). GSs-specific prognostic value, clinical characteristics, genes and pathways, immune infiltration, and chemotherapeutic drug sensitivity were evaluated. Three hub genes, KRT8, NGFR, and TCEA3, were screened and validated among GSs. These results suggested that the irlncRNAs had the promising potential for clinical immunotherapy of GBM.
Materials and methods
Acquisition and Processing of GBM Expression and Clinicopathological Data
RNA-seq transcriptome data of healthy samples were obtained from the GTEx database (11) (http://commonfund.nih.gov/GTEx/). The RNA-seq transcriptome data and clinicopathological data of the GBM samples were downloaded from the TCGA database (http://cancergenome.nih.gov/). Samples and patients with incomplete clinical information were excluded, and conformers are shown in Table S1. Two available matrices were merged, normalized with the limma package of R software, and obtained the differentially expressed (DE) genes. The input file is FPKM, and the output file is log (x + 1). The scRNA-seq data of human GBM samples, accession number GSE168004, were obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database. The cutoff criteria were set as | log2 fold change (FC) | > 0.5 and p < 0.05.
Identification of Immune-Related lncRNAs (irlncRNAs)
The immune genes (IGs) list was downloaded from the IMMPORT shared database (12) (https://www.immport.org/) and the Molecular Signatures Database v 7.0 (http://www.gsea-msigdb.org/gsea/index.jsp/). The correlation between genes was calculated to obtain irlncRNAs. Correlation coefficient >0.4 and p<0.001 were used as the threshold.
Establishment of the Immune-Related Risk Prognostic Model
Univariate and multivariate Cox regression analyses were performed to identify significant lncRNAs for construction of the prognostic signature. A risk score was calculated based on each patient’s lncRNAs expression level by the following formula: (N is the number of relative lncRNAs, Expi represents the expression value of each lncRNA, and βi is the regression coefficient of the multivariate Cox analysis for the target lncRNA). By setting the median value of the risk score as the cutoff value in the training set and the whole set, GBM patients were divided into high- and low-risk groups. Related files for constructing the immune-related risk prognostic model are displayed in Table S2.
Evaluation and Validation of a Risk Prognostic Model
The predictive ability of the prognostic model was evaluated by a series of analyses: Kaplan-Meier survival analysis, time-dependent ROC curve analyses, univariate Cox regression analysis, and multivariate Cox regression analysis for comparison of the survival between the high- and low-risk groups in the training, testing, and whole cohorts using the R packages survival and survivalROC. In addition, the signature derived from this study was compared with these three other signature ROC curves (13–15). We analyzed the ROC curve differences between prognostic models and clinicopathological features.
Construction of a CeRNA Network and a Protein–Protein Interaction (PPI) Network
The miRcode database (16) was performed to match differentially expressed and prognostically related irlncRNAs and miRNAs. Three databases, miRTarBase (17), miRDB (18), and TargetScan (19), were used to predict miRNA target genes. The interactions between miRNAs and lncRNAs or mRNAs were integrated to construct a CeRNA regulatory network. The mRNAs were enrolled in a PPI network through the STRING database (https://string-db.org/) with a confidence score > 0.7. Cytoscape (version 3.8.1) was used to visualize the CeRNA and PPI networks.
Functional and Pathway Enrichment Analyses
We used the “clusterProfiler” package to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of CeRNA network-related IGs to explore potential biological functions and pathways. The cutoff criterion was set at p < 0.05. Additionally, KEGG pathway analysis of KRT8 was performed using the Gene Set Enrichment Analysis (GSEA) software (www.gsea-msigdb.org). The cutoff criterion was set at p < 0.05.
Identification of GBM Subtypes in Risk Prognostic Model
Unsupervised consensus clustering was conducted to identify a novel immune classification of GBM based on the prognostic irlncRNAs using the ConsensusClusterPlus package (50 iterations, resample rate of 80%). The consensus cumulative distribution function (CDF), consensus matrix (CM), and consensus heatmap were performed to determine the optimal number of clusters.
Analysis of Clinical Characters and Molecular Differences in GBM Subtypes
Survival analysis and valuable clinical information (Table S3) were compared between the different subtypes. The Wilcoxon rank test was used to identify differentially expressed molecules among subtypes. The cutoff criteria were set as | log2FC | > 0.3 and p < 0.05.
Immune Microenvironment Exploration for GBM Subtypes Access scRNA-seq Data
The Seurat package was performed for quality control, statistical analysis, and exploration of the scRNA-seq data. The quality control standards were genes detected in >3 cells; cells with >50 total detected genes and cells with ≤5% of mitochondria-expressed genes were included. PCA was used to discriminate available dimensions with a p value < 0.05. Then, dimensionality reduction and cluster classification analysis were performed using a t-distributed stochastic neighbor embedding (tSNE) algorithm. The limma package was applied for differential expression analysis to identify the marker genes of each cluster with p value < 0.05 and | log2[fold change (FC)] | > 0.5. Based on marker gene populations, different cell clusters were annotated by the singleR package and then manually validated and corrected with the CellMarker database. The corresponding cell surface marker genes for the annotation of cell clusters are listed in Table S4.
Exploration of Candidate Small Molecule Agents
To evaluate the significance of this prognostic model in clinical treatment, the IC50 of common administrating chemotherapeutic agents in the GBM dataset TCGA project was calculated. The IC50 difference analysis was performed between the high-risk and low-risk groups using the Wilcoxon signed-rank test. Box plots were obtained using pRRophetic and ggplot2 to show the results.
Preparation for Human GBM Samples
GBM tissues and normal brain tissues were obtained from patients treated at First Affiliated Hospital of Nanchang University who provided informed consent. The study was approved by the hospital’s institutional ethics committee. GBM tissue was collected and immediately stored in an environment at −80°C.
Quantitative Real-Time RT-PCR (qRT-PCR) Analysis
Total RNA extracted from transfected cells was reverse-transcribed with RT reagent Kit gDNA Eraser (TaKaRa) and detected by SYBR-Green (TaKaRa). The PCR primers are listed in Table S5.
Western Blot Assays
Western blot (WB) assays were performed as described previously (20). The antibodies used are listed in Table S6.
Statistical Analysis
All analysis was carried out by R version 3.6.1 and corresponding packages. Kruskal–Wallis test was used to compare the divergence between multiple groups. Chi-square test or Fisher exact test was used for statistics on clinical information. A Bonferroni test was used to correct the p-value. Kaplan–Meier curves analysis was used to assess survival differences of the subtype. The correlation was determined by Pearson correlation analysis. p < 0.05 was regarded as statistically significant.
Results
Selection of DElncRNAs, DEimmune genes (DEIGs), and irlncRNAs in GBM
A filtering flow chart for the study is shown in Figure 1. The 1,520 DElncRNAs with 396 upregulated and 1,124 downregulated were identified between normal and GBM tissues. Analogously, we also identified 358 upregulated and 196 downregulated IGs. The corresponding heatmaps are displayed in Figure S2. Based on 396 upregulated lncRNAs and 554 DEIGs, 224 irlncRNAs were obtained by correlation analysis (Tables S7–S9).
Figure 1 Flow chart of study design.
Construction of irlncRNAs Model in GBM
The data of GBM patients were allocated randomly to the training and validation cohort. The 224 irlncRNAs were subjected to univariate Cox regression analysis (Table S10) followed by Lasso regression (Table S11 and Figures 2A–D) in the training set to obtain 17 PRirlncRNAs (p < 0.05; Table S12 and Figure 2E) and a risk score prognostic model constituted based on 6 key irlncRNAs. The risk score for each sample was calculated based on the expression levels of these six lncRNAs (Figure 2F). The coefficient of each gene was calculated by multivariate Cox regression analysis (Table 1).
Figure 2 Construction, evaluation, and comparison of a risk signature. (A) LASSO coefficient profiles of the 17 irlncRNAs in the training set. (B) A coefficient profile plot was generated against the log (lambda) sequence. Selection of the optimal parameter (lambda) in the LASSO model. (C, D) The AUC value and cutoff point obtained in the training set. (E) Forest plot of 17 irlncRNAs selected by univariate Cox regression analysis associated with GBM survival in the training set. (F) Forest plot of six irlncRNAs selected by multivariate Cox regression analysis associated with GBM survival and construction risk model. (G) Risk score and survival status analysis of irlncRNAs prognostic signature. (H) The expression pattern of irlncRNAs prognostic signature in the low- and high-risk groups. (I) Survival analysis of irlncRNAs prognostic signature. (J) ROC curve analysis within 1, 2, and 3 years. (K) Multivariate ROC curve analysis showing that the superior prognostic performance of the irlncRNAs prognostic model compared to other clinical indicators. (L) AUCs of the ROCs for our and the three other gene signatures.
Table 1 Coefficients based on a multivariate Cox regression analysis of the selected lncRNAs.
Evaluation and Validation of irlncRNAs Signature in GBM
The irlncRNAs signature is a robust prognostic tool for GBM. Risk curves and scatter plots showed the risk score and survival status of each GBM patient in the training (Figure 2G), testing (Figure S3A), and total sets (Figure S1A). The low-risk group had a lower risk coefficient and mortality than the high-risk group. The heatmap of the irlncRNAs signature in the training (Figure 2H; Table S13), testing (Figure S3B; Table S14), and total sets (Figure S1B) revealed that GBM with high prognostic scores expressed high-risk irlncRNAs (H19, AL162231.2, AC002456.1), whereas GBM with low prognostic scores expressed protective irlncRNAs (ST3GAL6-AS1, SOX21-AS1, AC006213.5). Based on the median risk score in the training set, GBM patients were divided into high- and low-risk cohorts. Survival curves indicated that patients in the low-risk group had a longer median OS compared with the high-risk group (Figure 2I); further examinations were performed in the test (Figure S3C) and whole sets (Figure S1C) by the same algorithmic cutoff in order to evaluate the accuracy of the prognostic signature. Both groups yielded similar results, suggesting that the prognostic signature was effective. In addition, the promising predictive value for the GBM special model in the training set was demonstrated by ROC curve analysis (Figure 2J 1-year AUC = 0.792, 2-year AUC = 0.922, 3-year AUC = 0.981), which validated the results of the model in the testing set (Figure S3D; 1-year AUC = 0.703, 2-year AUC = 0.657, 3-year AUC = 0.669) and the whole set (Figure S1D; 1-year AUC = 0.744, 2-year AUC = 0.756, 3-year AUC = 0.838). The multi-index ROC analysis revealed that the AUC of the prognostic model was significantly better than those of other clinicopathological indicators (Figure 2K) (such as age, gender, therapy, molecular typing, etc.). Compared with three existing lncRNA-related signatures (13–15) (Figure 2L), the excellent predictive viability of our model is further demonstrated. Together, these data illustrate the excellent identification of high-risk patients using our model.
IrlncRNAs Prognostic Model Is an Independent Prognostic Factor for GBM
Univariate and multivariate Cox regression analyses were performed to verify that the irlncRNAs model was an independent prognostic factor for GBM in the training set. The univariate Cox analysis revealed that gender, radiotherapy, MGMT status, and risk score were dramatically associated with the OS (Figure S4A), while the multivariate analysis revealed that gender, MGMT status, and risk score were identified as independent prognostic factors (Figure S4B).
Construction of the CeRNA and PPI Network and Functional Enrichment Analysis in GBM
Of the 224 irlncRNAs, 17 lncRNAs were associated with prognosis. Based on matching analysis of 17 PRirlncRNAs and 554 DEIGs, a total of 5 irlncRNAs and 16 miRNAs paired into 31 irlncRNAs–miRNA interactions, while 16 miRNAs and 27 DEIGs matched to form 35 miRNA–DEIGs interactions. Finally, 5 irlncRNAs, 16 miRNAs, and 27 DEIGs were used to construct lncRNA–miRNA–mRNA regulatory networks (Table S15 and Figure 3A).
Figure 3 GBM-specific CeRNA network, PPI network, and functional enrichment analysis (A) A total of 31 irlncRNAs–miRNA interactions and 35 miRNA–DEIGs interactions construct the lncRNA–miRNA–mRNA regulatory networks. (B) PPI network displayed the interactions of proteins translated from IGs in the CeRNA network. (C, D) GO enrichment analysis (E, F) KEGG pathway enrichment analysis.
Furthermore, the PPI network was constructed to identify the interactions of proteins translated from mRNAs in the CeRNA network (Figure 3B). We found that some genes with high combined score including TGFBR1-TGFBR2, JAG2-NOTCH2, ETS1-SP1, NRAS-PDGFRA, and BDNF-TRAF6 were mainly enriched in the “Human T-cell leukemia virus 1 infection,” “IL-17 signaling pathway,” “TGF-beta signaling pathway,” and “PD-L1 expression and PD-1 checkpoint pathway in cancer pathway.”
GO (Table S16) and KEGG (Table S17) pathway enrichment analyses demonstrated that the GBM-specific CeRNA network might be involved in the neoplastic process by regulating these biological functions and pathways. GO functional analysis showed that DEIGs involved in the CeRNA network were enriched in BPs, including regulation of vasculature development, response to oxidative stress, and positive regulation of epithelial cell proliferation. The enrichment of MF is mainly related to the membrane signal, and CC is protein binding (Figures 3C, D). CeRNA network-related IGs were significantly enriched in KEGG pathways, namely, MAPK signaling pathway, cytokine–cytokine receptor interaction, PI3K-Akt signaling pathway, and multiple cancers (Figures 3E–F).
Four Subtypes of GBM Were Identified and Correlated With Prognosis
Based on 17 PRirlncRNAs, Consensus Cluster Plus was utilized to identify the different subtypes (K = 2-9) among the risk model. According to the cumulative distribution function (CDF) curves, tracking pot, Delta area pot, and CM heatmap (Figures 4A–D), when k=4, the sample cluster was stable and robust. As a result, patients could be classified into four GSs (Table S18): A (n = 23, 27.4%), B (n = 24, 28.6%), C (n = 28, 33.3%), and D (n = 9, 10.7%). Kaplan-Meier survival analysis indicated that patients with GS-A showed the best OS compared to patients with cluster B, C, or D (p=0.007; Figure 4E).
Figure 4 Identification of potential GBM subtypes (A–C) Cumulative distribution function curve, delta area, and tracking plot of immune-related lncRNA in GBM. (D) Consensus clustering matrix for k = 4, which was the optimal cluster number in the TCGA training cohort. (E) Patients in the GBM subtype-A experienced a longer survival time.
Revelation of Clinical Characters, Molecular Differences, and Pathway Analysis for GBM Subtypes
Clinicopathological variables and molecular differences were compared among the four subtypes. Heatmap of 17 irlncRNAs illustrated the clinical features and molecular differences among the four subtypes (Figures 5A). The results revealed that GS-A patients are characterized by a high mutation rate of genes including IDH1, ATRX, and EGFR, a high rate of chemoradiotherapy, and a high rate of the low-risk group (Figures 5B–G).
Figure 5 Heatmap and clinicopathological features of four GBM subtypes (A) The heatmap and clinicopathological features of the 4 clusters based on the expression patterns of the 17 irlncRNAs in the training set. (B–F) Distribution ratio of IDH/ATRX/EGFR status and chemotherapy and radiotherapy in GBM subtypes. (G) Sankey diagram showing the prognosis of four GBM subtypes.
Subsequently, difference analysis identified 10 lncRNAs (Table S19) and 14 mRNAs (Table S20) among the four subtypes (Figures 6A, B). The results revealed that 6 of the 14 mRNAs were risk genes, and 4 (KRT8, NGFR, TCEA3, and PTTG1) of the risk genes were highly expressed in GBM compared with normal tissues. Thus, these four risk factors, as hub genes, may play an important role in the malignant behavior of GBM.
Figure 6 Molecular difference analysis of four GBM subtypes (A, B) Heatmaps of 10 differentially expressed irlncRNAs (A) and 14 mRNAs (B) between the 4 GBM subtypes. (C–E) GSEA showing that the functional pathways involved in RKT8 were mainly immune cell- and tumor-related signaling pathways. (F–H) GSVA enrichment analysis showing the activation states of biological pathways in GSs.
In addition, patients with GS-A patients are characterized by low expression of four high-risk lncRNA (H19, HOTAIRM1, AGAP2-AS1, AC002456.1) and one high-risk gene KRT8. GSEA showed that functional pathways involved in RTK8 were mainly immune cell and tumor-related signaling pathways, such as the T cell receptor, apoptosis, or JAK/STATA signaling pathway (Figures 6C–E).
GSVA enrichment analysis showed the activation states of biological pathways including the regulation of autophagy, the apoptosis, the Wnt signaling pathway, the NOTCH signaling pathway, the ERBB signaling pathway, the RIG like receptor signaling pathway, and the NOD-like receptor signaling pathway in GS-A (Figures 6F–H).
Exploration of Immune Microenvironment for GBM Subtypes Access scRNA-seq Data
Eight cell clusters with different annotations were identified by scRNA-seq data, revealing cellular heterogeneity in GBM tumors. A total of 4,210 cluster markers were identified from all 8 clusters by differential analysis (Table S21). Clusters 0 and 4, containing 398 cells, were annotated as GBM MSCs; clusters 1, 2, 3, 5, 6, and 7, containing 792 cells, were annotated as the astrocytes (Figures 7A–F and Table S22).
Figure 7 Estimation of tumor-infiltrating cells for GBM subtypes based on scRNA-seq data (A) After quality control of the 3,483 cells from the tumor cores of 4 human GBM samples, 1,190 cells were included in the analysis. (B) The variance diagram shows 13,859 corresponding genes throughout all cells from GBMs. The red dots represent highly variable genes, and the black dots represent nonvariable genes. The top 10 most variable genes are marked in the plot. (C) PCA identified the 15 PCs with an estimated p value < 0.05. (D) All eight clusters of cells in GBMs were annotated by singleR and CellMarker according to the composition of the marker genes. (E) The tSNE algorithm was applied for dimensionality reduction with the 20 PCs, and 8 cell clusters were successfully classified. (F) The differential analysis identified 4,210 marker genes. The top 20 marker genes of each cell cluster are displayed in the heatmap. A total of 68 genes are listed beside of the heatmap after omitting the same top marker genes among clusters. The colors from purple to yellow indicate the gene expression levels from low to high. (G–J) Expression profiles of the four risk genes in eight cell clusters.
Among the six risk genes of GBM subtypes, two genes, OLFM1 and TENM2, with low expression in GBM were excluded. Searching of the remaining four hub genes (GS-A: KRT8; GS-B: NGFR; GB-C: TCEA3; GB-D: PTTG1) in different GSs with clustering markers revealed that GBM may infiltrate immune cells. KRT8, belonging to Cluster 0, was annotated as GBM MSCs; NGFR, belonging to Clusters 0 and 4, was annotated as GBM MSCs; TCEA3, belonging to Clusters 0, 1, 2, 3, 4, and 5, was annotated as GBM MSCs and astrocyte; PTTG1, belonging to Clusters 2, 3, 4, 5, 6, and 7, was annotated as GBM MSCs and astrocyte (Figures 7G–J).
Screening of the Related Small Candidate Drugs With irlncRNAs Signature
An attempt was made to screen out chemotherapeutic agents that are sensitive to the high-risk group in the TCGA project of the GBM dataset. We found that the high-risk score correlated with a lower half inhibitory concentration (IC50) of chemotherapeutic drugs such as Gefitinib (p= 0.015) and Roscovitine (p= 0.035), whereas it correlated with the higher IC50 of axitinib (p=0.039) and thapsigargin (p=0.0041) (Figures S4A–D).
Validation of the Hub Genes in Clinical Tissues
The K-M survival curve from the TCGA database was performed to explore the potential role of the individual hub gene in OS. Three of the four hub genes showed significant predictions of poor OS (P < 0.05, Figures 8A–C). To further verify the expression level of hub genes in GBM samples, we generated RT-qPCR to calculate the mRNA levels of the three hub genes. As illustrated in Figures 8D–F, the expressions of KRT8, NGFR, and TCEA3 were significantly upregulated in GBM tissues compared with normal tissues. Subsequently, WB was used to evaluate the expression level of three proteins. As shown in Figures 8G–I, the expression levels of three proteins in GBM tissues were higher than those in normal brain tissues.
Figure 8 Validation of the hub genes in clinical tissues (A–C) Kaplan–Meier survival curves for patients of GBM with high and low gene expression in the TCGA dataset. (D–F) qRT-PCR of KRT8, NGFR, and TCEA3 in clinical human GBM samples and normal brain tissues. The expression levels were normalized to β-actin. ***P < 0.001. (G–I) The protein expression of KRT8, NGFR, and TCEA3 in clinical human GBM tissues and normal tissues was detected by WB.
Discussion
GBM cells form a complex tumor microenvironment that supports malignant tumor progression and immune escape (21). Novel immunotherapy within the tumor microenvironment has been uncovered that exerts antitumor immune response via targeting immunoregulatory cells or immunosuppressive factors (22). Accumulating evidence suggests that abnormal lncRNAs servers as new markers contribute to antitumor immunoreactivity (23, 24). It is of great significance to understand the tumor immune microenvironment driven by lncRNAs, to construct a clinical prognosis model, and to screen new markers for providing risk stratification and targets for immunotherapy.
In this study, 224 irlncRNAs were analyzed between tumor and normal tissues, 17 PRirlncRNAs were obtained by using the univariate Cox regression analysis, LASSO regression analysis was used to identify 6 key lncRNAs, and multivariate Cox regression analysis was applied to calculate coefficients and construct the risk model. We found that patients in the low-risk group had longer survival than those in the high-risk group. Subsequently, we established forest plots and ROC plots including age, sex, radiotherapy, chemotherapy, gene (IDH, MGMT, ATRX, and EGFR) mutation status, and risk scores. By plotting risk heatmap, risk curve, ROC curve, and survival curve, it was concluded that the risk model indeed had a good predictive effect. Meanwhile, we obtained similar results in the validation set.
The immune alterations driven by lncRNAs in GBM have also been preliminarily investigated (20, 25). Among the six key lncRNAs, H19, AL162231.2, and AC002456.1 were risk factors for the prognosis of GBM, while ST3GAL6-AS1, SOX21-AS1, and AC006213.5 were protective factors. LncRNA H19 as the first discovered classical regulator lncRNA is involved in the regulation of multiple cancers, including GBM (26, 27). H19 is overexpressed in glioma tissues, negatively correlates with patient survival, and promotes tumor growth by silencing relevant microRNAs (27, 28). H19 has a potential reference value for glioma remission and immunotherapy. ST3GAL6-AS1 and SOX21-AS1 as protective factors have been reported in cancers, lncRNA ST3GAL6-AS1 overexpression significantly reduces colorectal cancer cell tumorigenesis and metastasis (29), and lncRNA SOX21-AS1 significantly suppresses tumorigenesis in cervical cancer (30), oral cancer (31), and GBM (32). Relevant literature reports for AL162231.2, AC002456.1, and AC006213.5 are sparse.
Identifying the targets of lncRNAs is a key step in exploring their functions. An immune-related CeRNA network was constructed to predict lncRNAs targets, and a PPI network was used to evaluate the interactions of translated proteins from mRNAs in the CeRNA network. The CeRNA network enabled not only a deeper understanding of the communication between RNAs and a more comprehensive analysis of the complex gene interactions underlying carcinogenesis but also the identification of novel biomarkers. Among the prognostic biomarkers involved in the GBM-specific CeRNA and PPI network, the most significant difference gene was PLAU (encoding urokinase-type plasminogen activator; uPA), which was overexpressed, was the target of lncRNA H19, and was enriched in the KEGG pathway, namely, MicroRNAs in cancer, NF-kappa B signaling pathway, transcriptional misregulation in cancer, and proteoglycans in cancer pathways. Moreover, the protein pair with the highest combined score was TGFBR1–TGFBR2. PLAU is frequently upregulated in GBM (33, 34) and promotes cell invasion by PLAUR (PLAU receptor) binding and activation of extracellular proteases (35). TGFBR1 and TGFBR2 have been identified in GBM as a TGF‐β signaling upstream receptor (36), which has been well known to be a key regulator of migration phenotype in GBM cells (37). In our analysis, lncRNA H19 may exert biological activity by targeting miR-193a-3p to regulate gene PLAU expression. Therefore, it is meaningful to construct immune-related CeRNA and PPI networks in GBM to mine novel biomarkers, predict prognosis, and guide therapy.
In addition to identifying candidate biomarkers in GBM, GSs are also the key to improve personalized treatment (38). Based on 17 PRirlncRNAs, we can classify GBM patients into 4 GSs (A-D). Then, we assessed subtype-specific prognostic values, clinical characteristics, genes and pathways, immune infiltration, and chemotherapeutic drug sensitivity. Our results revealed that GS-A patients displayed the most favorable prognosis, which were characterized by a high mutation rate of genes including IDH1, ATRX and EGFR. Previously published reports indicated that IDH, ATRX and EGFR mutation status significantly influenced the prognosis of glioma patients (39). Such as, IDH mutations are frequent in infiltrating astrocytomas (grades II and III) and secondary GBMs (1). Primary GBMs typically lack IDH mutations and demonstrate EGFR, PDGFRA, TP53, PTEN, NF1, and TERT promoter mutations (40). These classical biomarkers have been integrated into multiple classification schemes and applied to an accurate clinical decision-making process. We observed that GBM with GS-A were characterized by four high-risk lncRNAs (H19, HOTAIRM1, AGAP2-AS1, and AC002456.1) and one high-risk mRNA KRT8 with a low expression level. Among these lncRNAs and mRNAs, lncRNA HOXA transcript antisense RNA myeloid-specific 1 (HOTAIRM1) participates in the reprogramming of chromatin organization and proliferation and metastasis of cancer (32), which has been found to be highly expressed in a variety of tumors including GBM (41). LncRNA AGAP2 antisense RNA 1 (AGAP2-AS1), transcribed from a gene located at 12q14.1, a novel cancer-related lncRNA, was dysregulated in cancers (42). In GBM, lncRNA HOTAIRM1 (43, 44) and AGAP2-AS1 (45, 46), as oncogenic factors, promoted tumorigenesis, predicted a poor clinical outcome, and were potential biomarkers and therapeutic targets. Keratin 8 (KRT8), a major component of the intermediate filament cytoskeleton, promotes tumor progression and metastasis of various cancers (47–49). In our analysis, GS-A was positively correlated with autophagy, the apoptosis, the Wnt signaling pathway, the NOTCH signaling pathway, the ERBB signaling pathway, the RIG like receptor signaling pathway, and the NOD-like receptor signaling pathway. KRT8 may exert its biological activity through regulating the T cell receptor, apoptosis, or the JAK/STATA signaling pathway.
The efficacy of immunotherapy strongly depends on intertumoral tumor-infiltrating immune cells (50). Combining the risk gene of GSs and scRNA-seq data reveals poor prognosis GBM tumor-infiltrating immunoreactive cells. Park et al. demonstrated that high expression on macrophage signatures of GBM patients predicted suboptimal survival (51), which was consistent with our analysis. We observed that the major tumor-infiltrating immune cells in GBM with poor prognosis were MSCs and astrocyte. Radfar et al. demonstrated that nonspecific activation of CD4+ T cells dramatically enhanced the cytotoxicity of four chemotherapeutic agents including TMZ, paclitaxel (Pax), Carbo, and 5-FU in cancers (52). Patients with more infiltrated CD8+ T cells had a better response to pembrolizumab treatment than those with less infiltrated cells (53). Our model suggested that high risk was associated with sensitivity to chemotherapeutics such as gefitinib and roscovitine, and GS-A patients in low risk were more sensitive to axitinib and thapsigargin.
Our survival analysis and in vitro study showed that three of the four hub genes showed significant predictions of poor OS, and the mRNA and protein levels of KRT8, NGFR, and TCEA3 were significantly upregulated in GBM tissues compared with normal tissues. The above results indicated that these proteins encoded by the hub genes may play a feasible oncogenic role in GBM.
Conclusion
In this study, based on 17 PRirlncRNAs, we not only constructed a six-key irlncRNAs prognostic signature but also identified four subtypes of GBM, which had a potential prognostic value. In GBM, lncRNA H19 may exert biological activity by targeting miR-193a-3p to regulate gene PLAU expression; KRT8, NGFR, and TCEA3 may stimulate novel strategies for immunotherapy of GBM patients. Interestingly, KRT8 may exert its biological activity through regulating the T cell receptor, apoptosis, or the JAK/STATA signaling pathway.
Highlights
1. Transcriptome and clinical information from 168 GBM samples was employed to screen 17 immune related lncRNAs (irlncRNAs) associated with prognosis.
2. 17 PRirlncRNAs were screen to construct a signature of 6 key irlncRNAs, which showing a good predictive effect, and similar results in the validation set.
3. GBM-specific immune CeRNA and PPI networks were constructed to predict lncRNAs targets and evaluate the interactions and functions of immune mRNAs translated proteins based on 17 PRirlncRNAs.
4. Four GBM subtypes (A–D) were identified based on 17 PRirlncRNAs, and we evaluated subtype-specific prognostic values, clinical characteristics, genes and pathways, immune infiltration, and chemotherapeutics efficacy.
5. Construction of the lncRNAs risk model and identification of GBM subtypes under immune environment, suggesting that KRT8, NGFR, TCEA3, and irlncRNAs had promising potential for clinical immunotherapy of GBM.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
Consent for participation for all patients was obtained through The Genotype-Tissue Expression (GTEx) Project, The Cancer Genome Atlas (TCGA) Project, and the Gene Expression Omnibus (GEO) database. Paraffin-embedded GBM tissues and normal brain tissues were obtained from patients who provided informed consent under an Institutional Ethics Committee-approved study from the First Affiliated Hospital of Nanchang University.
Author Contributions
All authors contributed to the analysis of data in this study. Conception and design: GZ, FQ and PW. Acquisition, analysis and interpretation of data: WY, WH and FW; Writing, review, and/or revision of the manuscript: WY; Administrative, technical, or material support: YM and PW. Study supervision: YM and WC. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the Project of Nanchang Science and Technology Support Plan of Jiangxi Province, China (HONGKO Zi [2021] 129) and the National Institute of Health Grant (5R01AR067319-04. 09/2015-07/2020).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
The reviewers are grateful for their helpful comments on this article.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.706936/full#supplementary-material
Abbreviations
GBM, glioblastoma multiforme; lncRNAs, long noncoding RNAs; irlncRNAs, immune-related lncRNAs; GSs, GBM subtypes; PRirlncRNAs, prognostically related irlncRNAs; GTEx, genotype-tissue expression; TCGA, The Cancer Genome Atlas; CGGA, Chinese Glioma Genome Atlas; GEO, Gene Expression Omnibus; OS, overall survival; IGs, immunity genes; NEAT1, nuclear-enriched abundant transcript 1; TIM, tumor immune microenvironment; DE, differentially expressed; IG, immune genes; PPI, protein–protein interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; ROC, receiver operating characteristic; TMA, temozolomide; FC, fold change; CDF, cumulative distribution function; CM, consensus matrix; BP, biological processes; MF, molecular function; CC, cellular component; IDH, isocitrate dehydrogenase; MGMT, O6-methylguanine DNA methyltransferase; EGFR, epidermal growth factor receptor; HOTAIRM1, HOXA transcript antisense RNA myeloid-specific 1; AGAP2-AS1, AGAP2 antisense RNA 1; KRT8, keratin 8; NGFR, nerve-growth-factor receptor; TCEA3, transcription elongation factor A3; PTTG1, pituitary tumor transforming gene 1.
References
1. Aldape K, Zadeh G, Mansouri S, Reifenberger G, von Deimling A. Glioblastoma: Pathology, Molecular Mechanisms and Markers. Acta Neuropathol (2015) 129:829–48. doi: 10.1007/s00401-015-1432-1
2. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, et al. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: A Summary. Acta Neuropathol (2016) 131:803–20. doi: 10.1007/s00401-016-1545-1
3. Bhan A, Soleimani M, Mandal SS. Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res (2017) 77:3965–81. doi: 10.1158/0008-5472.CAN-16-2634
4. Martens-Uzunova ES, Bottcher R, Croce CM, Jenster G, Visakorpi T, Calin GA. Long Noncoding RNA in Prostate, Bladder, and Kidney Cancer. Eur Urol (2014) 65:1140–51. doi: 10.1016/j.eururo.2013.12.003
5. Statello L, Guo CJ, Chen LL, Huarte M. Gene Regulation by Long non-Coding RNAs and its Biological Functions. Nat Rev Mol Cell Biol (2021) 22:96–118. doi: 10.1038/s41580-020-00315-9
6. Mathy NW, Chen XM. Long non-Coding RNAs (lncRNAs) and Their Transcriptional Control of Inflammatory Responses. J Biol Chem (2017) 292:12375–82. doi: 10.1074/jbc.R116.760884
7. Hirose T, Virnicchi G, Tanigawa A, Naganuma T, Li R, Kimura H, et al. NEAT1 Long Noncoding RNA Regulates Transcription via Protein Sequestration Within Subnuclear Bodies. Mol Biol Cell (2014) 25:169–83. doi: 10.1091/mbc.e13-09-0558
8. Hu G, Gong AY, Wang Y, Ma S, Chen X, Chen J, et al. LincRNA-Cox2 Promotes Late Inflammatory Gene Transcription in Macrophages Through Modulating SWI/SNF-Mediated Chromatin Remodeling. J Immunol (2016) 196:2799–808. doi: 10.4049/jimmunol.1502146
9. Zhao L, Liu Y, Zhang J, Liu Y, Qi Q. LncRNA SNHG14/miR-5590-3p/ZEB1 Positive Feedback Loop Promoted Diffuse Large B Cell Lymphoma Progression and Immune Evasion Through Regulating PD-1/PD-L1 Checkpoint. Cell Death Dis (2019) 10:731. doi: 10.1038/s41419-019-1886-5
10. Hu Q, Ye Y, Chan LC, Li Y, Liang K, Lin A, et al. Oncogenic lncRNA Downregulates Cancer Cell Antigen Presentation and Intrinsic Tumor Suppression. Nat Immunol (2019) 20:835–51. doi: 10.1038/s41590-019-0400-7
11. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: A Web Server for Cancer and Normal Gene Expression Profiling and Interactive Analyses. Nucleic Acids Res (2017) 45:W98–W102. doi: 10.1093/nar/gkx247
12. Shen S, Wang G, Zhang R, Zhao Y, Yu H, Wei Y, et al. Development and Validation of an Immune Gene-Set Based Prognostic Signature in Ovarian Cancer. EBioMedicine (2019) 40:318–26. doi: 10.1016/j.ebiom.2018.12.054
13. Wang W, Zhao Z, Yang F, Wang H, Wu F, Liang T, et al. An Immune-Related lncRNA Signature for Patients With Anaplastic Gliomas. J Neurooncol (2018) 136:263–71. doi: 10.1007/s11060-017-2667-6
14. Wang X, Gao M, Ye J, Jiang Q, Yang Q, Zhang C, et al. An Immune Gene-Related Five-lncRNA Signature for to Predict Glioma Prognosis. Front Genet (2020) 11:612037. doi: 10.3389/fgene.2020.612037
15. Xia P, Li Q, Wu G, Huang Y. An Immune-Related lncRNA Signature to Predict Survival In Glioma Patients. Cell Mol Neurobiol (2021) 41:365–75. doi: 10.1007/s10571-020-00857-8
16. Herbst RS, Heymach JV, Lippman SM. Lung Cancer. N Engl J Med (2008) 359:1367–80. doi: 10.1056/NEJMra0802714
17. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2019. CA Cancer J Clin (2019) 69:7–34. doi: 10.3322/caac.21551
18. Lin JJ, Cardarella S, Lydon CA, Dahlberg SE, Jackman DM, Janne PA, et al. Five-Year Survival in EGFR-Mutant Metastatic Lung Adenocarcinoma Treated With EGFR-TKIs. J Thorac Oncol (2016) 11:556–65. doi: 10.1016/j.jtho.2015.12.103
19. Slack FJ, Chinnaiyan AM. The Role of Non-Coding RNAs in Oncology. Cell (2019) 179:1033–55. doi: 10.1016/j.cell.2019.10.017
20. Wu P, Cai J, Chen Q, Han B, Meng X, Li Y, et al. Lnc-TALC Promotes O(6)-Methylguanine-DNA Methyltransferase Expression via Regulating the C-Met Pathway by Competitively Binding With miR-20b-3p. Nat Commun (2019) 10:2045. doi: 10.1038/s41467-019-10025-2
21. Hanahan D, Coussens LM. Accessories to the Crime: Functions of Cells Recruited to the Tumor Microenvironment. Cancer Cell (2012) 21:309–22. doi: 10.1016/j.ccr.2012.02.022
22. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the Tumor Microenvironment: Removing Obstruction to Anticancer Immune Responses and Immunotherapy. Ann Oncol (2016) 27:1482–92. doi: 10.1093/annonc/mdw168
23. Li Y, Jiang T, Zhou W, Li J, Li X, Wang Q, et al. Pan-Cancer Characterization of Immune-Related lncRNAs Identifies Potential Oncogenic Biomarkers. Nat Commun (2020) 11:1000. doi: 10.1038/s41467-020-14802-2
24. Atianand MK, Caffrey DR, Fitzgerald KA. Immunobiology of Long Noncoding RNAs. Annu Rev Immunol (2017) 35:177–98. doi: 10.1146/annurev-immunol-041015-055459
25. Peng Z, Liu C, Wu M. New Insights Into Long Noncoding RNAs and Their Roles in Glioma. Mol Cancer (2018) 17:61. doi: 10.1186/s12943-018-0812-2
26. Brannan CI, Dees EC, Ingram RS, Tilghman SM. The Product of the H19 Gene may Function as an RNA. Mol Cell Biol (1990) 10:28–36. doi: 10.1128/MCB.10.1.28
27. Jia P, Cai H, Liu X, Chen J, Ma J, Wang P, et al. Long non-Coding RNA H19 Regulates Glioma Angiogenesis and the Biological Behavior of Glioma-Associated Endothelial Cells by Inhibiting microRNA-29a. Cancer Lett (2016) 381:359–69. doi: 10.1016/j.canlet.2016.08.009
28. Deng Y, Zhou L, Yao J, Liu Y, Zheng Y, Yang S, et al. Associations of lncRNA H19 Polymorphisms at MicroRNA Binding Sites With Glioma Susceptibility and Prognosis. Mol Ther Nucleic Acids (2020) 20:86–96. doi: 10.1016/j.omtn.2020.02.003
29. Hu J, Shan Y, Ma J, Pan Y, Zhou H, Jiang L, et al. LncRNA ST3Gal6-AS1/ST3Gal6 Axis Mediates Colorectal Cancer Progression by Regulating Alpha-2,3 Sialylation via PI3K/Akt Signaling. Int J Cancer (2019) 145:450–60. doi: 10.1002/ijc.32103
30. Wang R, Li Y, Du P, Zhang X, Li X, Cheng G. Hypomethylation of the lncRNA SOX21-AS1 has Clinical Prognostic Value in Cervical Cancer. Life Sci (2019) 233:116708. doi: 10.1016/j.lfs.2019.116708
31. Yang CM, Wang TH, Chen HC, Li SC, Lee MC, Liou HH, et al. Aberrant DNA Hypermethylation-Silenced SOX21-AS1 Gene Expression and its Clinical Importance in Oral Cancer. Clin Epigenet (2016) 8:129. doi: 10.1186/s13148-016-0291-5
32. Wang Z, Ji X, Gao L, Guo X, Lian W, Deng K, et al. Comprehensive In Silico Analysis of a Novel Serum Exosome-Derived Competitive Endogenous RNA Network for Constructing a Prognostic Model for Glioblastoma. Front Oncol (2021) 11:553594. doi: 10.3389/fonc.2021.553594
33. Bernhart E, Damm S, Wintersperger A, DeVaney T, Zimmer A, Raynham T, et al. Protein Kinase D2 Regulates Migration and Invasion of U87MG Glioblastoma Cells In Vitro. Exp Cell Res (2013) 319:2037–48. doi: 10.1016/j.yexcr.2013.03.029
34. Dikshit B, Irshad K, Madan E, Aggarwal N, Sarkar C, Chandra PS, et al. FAT1 Acts as an Upstream Regulator of Oncogenic and Inflammatory Pathways, via PDCD4, in Glioma Cells. Oncogene (2013) 32:3798–808. doi: 10.1038/onc.2012.393
35. Blasi F, Carmeliet P. uPAR: A Versatile Signalling Orchestrator. Nat Rev Mol Cell Biol (2002) 3:932–43. doi: 10.1038/nrm977
36. Vinchure OS, Sharma V, Tabasum S, Ghosh S, Singh RP, Sarkar C, et al. Polycomb Complex Mediated Epigenetic Reprogramming Alters TGF-Beta Signaling via a Novel EZH2/miR-490/TGIF2 Axis Thereby Inducing Migration and EMT Potential in Glioblastomas. Int J Cancer (2019) 145:1254–69. doi: 10.1002/ijc.32360
37. Paw I, Carpenter RC, Watabe K, Debinski W, Lo HW. Mechanisms Regulating Glioma Invasion. Cancer Lett (2015) 362:1–7. doi: 10.1016/j.canlet.2015.03.015
38. Berger AC, Korkut A, Kanchi RS, Hegde AM, Lenoir W, Liu W, et al. A Comprehensive Pan-Cancer Molecular Study of Gynecologic and Breast Cancers. Cancer Cell (2018) 33:690–705 e699. doi: 10.1016/j.ccell.2018.03.014
39. van Tellingen O, Yetkin-Arik B, de Gooijer MC, Wesseling P, Wurdinger T, de Vries HE. Overcoming the Blood-Brain Tumor Barrier for Effective Glioblastoma Treatment. Drug Resist Update (2015) 19:1–12. doi: 10.1016/j.drup.2015.02.002
40. Appin CL, Brat DJ. Biomarker-Driven Diagnosis of Diffuse Gliomas. Mol Aspects Med (2015) 45:87–96. doi: 10.1016/j.mam.2015.05.002
41. Pastori C, Kapranov P, Penas C, Peschansky V, Volmar CH, Sarkaria JN, et al. The Bromodomain Protein BRD4 Controls HOTAIR, a Long Noncoding RNA Essential for Glioblastoma Proliferation. Proc Natl Acad Sci USA (2015) 112:8326–31. doi: 10.1073/pnas.1424220112
42. Hellner KA. [Night Blindness and Retinopathy in Systemic Diseases]. Verh Dtsch Ges Inn Med (1976) 82 Pt 1:580–3.
43. Li Q, Dong C, Cui J, Hong X. Over-Expressed lncRNA HOTAIRM1 Promotes Tumor Growth and Invasion Through Up-Regulating HOXA1 and Sequestering G9a/EZH2/Dnmts Away From the HOXA1 Gene in Glioblastoma Multiforme. J Exp Clin Cancer Res (2018) 37:265. doi: 10.1186/s13046-018-0941-x
44. Xie P, Li X, Chen R, Chen Y, Liu D, Liu W, et al. Upregulation of HOTAIRM1 Increases Migration and Invasion by Glioblastoma Cells. Aging (Albany NY) (2020) 13:2348–64. doi: 10.18632/aging.202263
45. Luo W, Li X, Song Z, Zhu X, Zhao S. Long non-Coding RNA AGAP2-AS1 Exerts Oncogenic Properties in Glioblastoma by Epigenetically Silencing TFPI2 Through EZH2 and LSD1. Aging (Albany NY) (2019) 11:3811–23. doi: 10.18632/aging.102018
46. Luo X, Tu T, Zhong Y, Xu S, Chen X, Chen L, et al. AGAP2-AS1 May Promote the Occurrence and Development of Glioblastoma by Sponging miR-9-5p: Evidence From a ceRNA Network. Front Oncol (2021) 11:607989. doi: 10.3389/fonc.2021.607989
47. Ramanand SG, Chen Y, Yuan J, Daescu K, Lambros MB, Houlahan KE. The Landscape of RNA Polymerase II-Associated Chromatin Interactions in Prostate Cancer. J Clin Invest (2020) 130:3987–4005. doi: 10.1172/JCI134260
48. Fang J, Wang H, Liu Y, Ding F, Ni Y, Shao S, et al. High KRT8 Expression Promotes Tumor Progression and Metastasis of Gastric Cancer. Cancer Sci (2017) 108:178–86. doi: 10.1111/cas.13120
49. Xie L, Dang Y, Guo J, Sun X, Xie T, Zhang L, et al. High KRT8 Expression Independently Predicts Poor Prognosis for Lung Adenocarcinoma Patients. Genes (Basel) (2019) 10. doi: 10.3390/genes10010036
50. McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, et al. Clonal Neoantigens Elicit T Cell Immunoreactivity and Sensitivity to Immune Checkpoint Blockade. Science (2016) 351:1463–9. doi: 10.1126/science.aaf1490
51. Iglesia MD, Parker JS, Hoadley KA, Serody JS, Perou CM, Vincent BG. Genomic Analysis of Immune Cell Infiltrates Across 11 Tumor Types. J Natl Cancer Inst (2016) 108. doi: 10.1093/jnci/djw144
52. Radfar S, Wang Y, Khong HT. Activated CD4+ T Cells Dramatically Enhance Chemotherapeutic Tumor Responses In Vitro and In Vivo. J Immunol (2009) 183:6800–7. doi: 10.4049/jimmunol.0901747
Keywords: glioblastoma, immune-related lncRNAs, biomarker, prognostic signature, immune infiltration
Citation: Yu W, Ma Y, Hou W, Wang F, Cheng W, Qiu F, Wu P and Zhang G (2021) Identification of Immune-Related lncRNA Prognostic Signature and Molecular Subtypes for Glioblastoma. Front. Immunol. 12:706936. doi: 10.3389/fimmu.2021.706936
Received: 08 May 2021; Accepted: 25 October 2021;
Published: 25 November 2021.
Edited by:
Mingzhu Yin, Central South University, ChinaReviewed by:
Zongqiang Huang, First Affiliated Hospital of Zhengzhou University, ChinaHongzhe Guo, Harbin Institute of Technology, China
Tao Zhennan, Nanjing Drum Tower Hospital, China
Copyright © 2021 Yu, Ma, Hou, Wang, Cheng, Qiu, Wu and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guohua Zhang, zhgh71@163.com; Pengfei Wu, wupengfei1992@126.com; Feng Qiu, lukeqiubmu@163.com
†These authors have contributed equally to this work