Molecular Subtypes Based on the Stemness Index Predict Prognosis in Glioma Patients

Glioma is the common histological subtype of malignancy in the central nervous system, with high morbidity and mortality. Glioma cancer stem cells (CSCs) play essential roles in tumor recurrence and treatment resistance. Thus, exploring the stem cell-related genes and subtypes in glioma is important. In this study, we collected the RNA-sequencing (RNA-seq) data and clinical information of glioma patients from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases. With the differentially expressed genes (DEGs) and weighted gene correlation network analysis (WGCNA), we identified 86 mRNA expression-based stemness index (mRNAsi)-related genes in 583 samples from TCGA RNA-seq dataset. Furthermore, these samples from TCGA database could be divided into two significantly different subtypes with different prognoses based on the mRNAsi corresponding gene, which could also be validated in the CGGA database. The clinical characteristics and immune cell infiltrate distribution of the two stemness subtypes are different. Then, functional enrichment analyses were performed to identify the different gene ontology (GO) terms and pathways in the two different subtypes. Moreover, we constructed a stemness subtype-related risk score model and nomogram to predict the prognosis of glioma patients. Finally, we selected one gene (ETV2) from the risk score model for experimental validation. The results showed that ETV2 can contribute to the invasion, migration, and epithelial-mesenchymal transition (EMT) process of glioma. In conclusion, we identified two distinct molecular subtypes and potential therapeutic targets of glioma, which could provide new insights for the development of precision diagnosis and prognostic prediction for glioma patients.

Introduction

Gliomas are the most common primary malignant tumor of the central nervous system (Ostrom et al., 2014). Therapeutic strategies, including surgery, chemotherapy, and radiotherapy, have been widely applied, but the overall outcome of glioma patients is still unsatisfactory (Bush et al., 2017; Nabors et al., 2017). The heterogeneity of glioma is an account of the poor prognosis (Klughammer et al., 2018). Thus, exploring the molecular mechanism of glioma may facilitate the identification of prognostic biomarkers and potential targets for the treatment of glioma.

Cancer stem cells (CSCs) are a subset of cancer cells with characteristics such as the ability to self-renew long-term, differentiate into defined progenies, and sustain tumor growth (Vlashi and Pajonk, 2015; Clarke, 2019). CSCs contribute to glioma recurrence, radioresistance, and chemoresistance through multiple molecular mechanisms (Prieto-Vila et al., 2017; Schulz et al., 2019). Therefore, a further understanding of the biological behavior of glioma stem cells may facilitate to changes in the treatment dilemma of glioma. Recently, the mRNA expression-based stemness index (mRNAsi), which represents the transcriptomic stemness expression, has been applied to assess CSC characteristics (Malta et al., 2018). In some cancers, such as bladder, lung, breast, or endometrial carcinoma, it has been reported that mRNAsi is a credible marker and is associated with tumor stage (Pan et al., 2019; Liu et al., 2020; Pei et al., 2020). However, its roles in glioma are rare known.

In this study, we collected the RNA-sequencing (RNA-seq) data and clinical information of glioma patients from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases. We found that glioma patients could be divided into two significant stemness subtypes (S1 and S2 groups) based on mRNAsi-related genes. We identified that the clinical characteristics, such as age, IDH status, and WHO grades, were different in the S1 and S2 groups. The tumor microenvironments in the two groups were also different. Furthermore, based on differentially expressed genes (DEGs) between the two stemness subtypes, a prognostic prediction model was constructed and could effectively divide patients into different prognoses in both TCGA and CGGA datasets. Finally, we selected one gene (ETV2) from the risk score model for further experimental validation. The results showed that ETV2 was more highly expressed in glioma and contributed to the invasion, migration, and epithelial-mesenchymal transition (EMT) process of glioma. Thus, our study provides novel molecular subtypes based on the stemness index to predict prognosis in glioma patients who may promote the clinical diagnosis and treatment of glioma.

Materials and Methods

Data Preparation

Figure 1 shows the workflow of data analysis. The RNA-seq data and corresponding clinical data of glioma (lower-grade glioma, LGGs and glioblastoma, GBM) and normal samples in TCGA database were downloaded from the UCSC Xena database¹ (dataset ID: TCGA.GBMLGG.sampleMap/HiSeqV2). In addition, another glioma dataset (RNA-seq and clinical data) was downloaded from the CGGA² (dataset ID: mRNAseq_693). The mRNAsi indices of glioma in TCGA were obtained from a previous study (Malta et al., 2018). The tumor purity of glioma is calculated based on the Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm, which can predict the level of infiltrating stromal and immune cells and then infer the tumor purity of tumor tissue (Yoshihara et al., 2013).

FIGURE 1

Figure 1. Flowchart presenting the process of analysis.

Differentially Expressed Genes Analysis

The “limma” R package was utilized to perform the DEG analysis between glioma and normal samples in TCGA database. |Log₂(Fold change)| > 1 and false discovery rate (FDR) < 0.05 were considered as the cutoff criteria. The volcano plot was drawn to show the DEGs.

Weighted Gene Correlation Network Analysis

The selected DEGs were used in weighted gene correlation network analysis (WGCNA) by the WGCNA R package (Langfelder and Horvath, 2008). After filtering the outliers in RNA-seq data, a Pearson correlation matrix was constructed for paired genes. Then, we established a weighted adjacency matrix by the power function amn = |cmn|β, as a previous study described (Yu et al., 2012). A proper β value was selected to increase the matrix similarity and establish a co-expression network. Next, the adjacency matrix was converted into a topological overlap matrix (TOM) to measure gene connectivity in the network. Based on TOM-based dissimilarity measurements, average linkage hierarchical clustering was performed with a minimum gene dendrogram size over 30. Finally, their dissimilarity was calculated, and module dendrograms were constructed for further analysis.

Gene significance (GS), representing the correlation between genes and sample traits, was calculated for each module. In the principal component analysis of each module, module eigengenes (MEs) were considered as the first principal component of a clustered module representing the gene expression profiles. Module membership (MM) was defined as the correlation between the module genes and gene expression profiles. In our study, mRNAsi and epigenetically regulated mRNAsi (EGER-mRNAsi) were the selected clinical phenotypes for further analysis. When GS and MM were highly correlated (more than 0.7), the module’s own genes were considered as significantly correlated with clinical traits.

Consensus Clustering

After finding the genes highly correlated with mRNAsi and EGER-mRNAsi, we performed consensus clustering to divide patients into different stemness subtypes based on these genes. The R package “ConsensusClusterPlus” was adopted to perform the consensus clustering (Wilkerson and Hayes, 2010). The optimal number of subgroups was determined by the cumulative distribution function (CDF) and consensus matrices. In addition, the CGGA dataset was used to validate this clustering.

Protein-Protein Interaction Analysis

To explore the protein-protein interaction (PPI) network of the selected genes, they were imported into the STRING database,³ which is a web tool used to explore the interactions between multiple proteins.

Comparison of Stemness Subtypes

To investigate the difference between the two stemness subtypes in TCGA and CGGA datasets, we used two-sample t-tests to compare different clinical variables such as age, IDH status, and WHO grades in the subtypes. Moreover, stratified survival analysis in different WHO grades was also used to evaluate the prognostic predictive value of the stemness subtypes.

Profiling of Immune Infiltrates in the Two Stemness Subtypes

Becht et al. (2016) designed the microenvironment cell populations-counter (MCP-counter) method, which can robustly quantify the absolute abundance of eight immune (T cells, CD8 T cells, cytotoxic lymphocytes, NK cells, B lineage, monocytic lineage, myeloid dendritic cells, and neutrophils) and two stromal cell (endothelial cells and fibroblasts) populations in heterogeneous tissues from transcriptomic data. The MCPs were analyzed by the “MCPcounter” R package. Based on this package, the proportions of eight tumor-infiltrating immune cells (T cells, CD8 T cells, cytotoxic lymphocytes, NK cells, B lineage, monocytic lineage, myeloid dendritic cells, and neutrophils) and two non-immune cells (endothelial cells and fibroblasts) were calculated based on the normalized RNA-seq data in TCGA and CGGA databases.

Functional Enrichment Analysis

First, DEG analysis between two stemness subtypes was performed. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID) to explore the different mechanisms and pathways between the two subtypes (Huang et al., 2009a,b).

Gene set enrichment analysis (GSEA) was performed between the two stemness subtypes with GSEA software. The reference gene set (c2.cp.kegg.v7.1.symbols.gmt) was acquired from the MSigDB database.⁴ Only enriched KEGG pathways with a p < 0.05 and FDR < 0.25 were considered statistically significant.

Construction and Validation of a Stemness Subtype-Related Prognostic Signature

First, DEG analysis between two stemness subtypes was performed in TCGA dataset as previously described. Then, univariate Cox hazard analysis was performed with p < 0.05 as a threshold parameter for the DEGs between the two stemness subtypes. By applying the “glmnet” R package, the least absolute shrinkage and selection operator (LASSO) algorithm was used to construct the stemness subtype-related prognostic signature (Goeman, 2010). This signature was created utilizing Cox regression coefficients to multiply the expression values of the select genes. According to the median of risk score, the patients were divided into high and low-risk groups. In addition, the time-dependent ROC curve and Kaplan-Meier survival curve analyses were used to evaluate the accuracy of the signature. The signature was also validated with CGGA dataset.

Nomogram Construction

To evaluate the predictive value of the prognostic signature, univariate and multivariate Cox regression analyses were performed together with the clinical information (grade, age, sex, IDH1, and 1p/19q status). Then, a nomogram was constructed to predict the survival probability by using the “rms” R package (Yin et al., 2020). Calibrations were used to evaluate the accuracy of the nomogram.

Cell Culture, Real-Time Quantitative PCR

HEB, SHG44, and A172 glioma cells were provided by Xiangya Medical School of Central South University, Changsha, China. HEB, A172, and SHG44 cells were cultured in DMEM high glucose medium (Gibco/Thermo Fisher Scientific, Inc.) with 10% fetal bovine serum at 37°C, 5%CO₂. SiRNAs were transfected with Lipofectamine 2000 (Thermo Fisher Scientific) 48 h before analysis. The siRNAs against the ETV2 gene were synthesized by RiboBio Corporation (Product number: siG000002116A-1-5, Guangzhou, China). Total RNA from HEB, SHG44, and A172 cells was extracted by the Trizol lysis method. cDNA synthesis was performed according to the Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA). The RNA levels of ETV2 were detected by using real-time quantitative PCR (qRT-PCR) according to the manufacturer’s protocol. The expression of ETV2 and GAPDH was analyzed by the 2^-ΔΔCt method. The primers were obtained from Sangon (Shanghai, China) and the sequences were designed as follows: for ETV2, the forward primer was 5'-CTGGAAAGGTACAAGCTCATCC-3' and the reverse primer was 5'-AACTTCTGGGTGCAGTAACGC-3'. For GAPDH, the forward primer was 5'-CATTGACCTCAACTACATGGTT-3' and the reverse primer was 5'-CCATTGATGACAAGCTTCCC-3'.

Wound Healing and Transwell assay

Wound healing and transwell assays were performed as previously described (Jia et al., 2018).

Western Blots

A172 and SHG44 cells were lysed with RIPA buffer for half an hour at 4°C. The supernatant was collected and boiled at 95°C for 5–8 min in SDS loading buffer. Then, they were subjected to electrophoresis in 10% SDS-polyacrylamide gels and transferred to the polyvinylidene difluoride membranes. The membranes were blocked with 5% non-fat milk in phosphate-buffered saline (PBS) for 1 h before being incubated with the primary antibody at 4°C overnight. The primary antibodies for western blotting used in this study were GAPDH, ETV2 (ab181847, Abcam), N-cadherin (22018-1-AP, proteintech), and vimentin (10366-1-AP, Proteintech). Then the cells were washed three to four times with 0.1% PBST and incubated with horseradish peroxidase (HRP)-conjugated secondary antibody (1:10,000) for 1 h at room temperature. The membranes were washed in 0.1% PBST four times before exposure. Chemiluminescent HRP substrate was purchased from Millipore (Catalog: WBKLS0500). Images were acquired in a Bio-Rad Universal Hood II machine with Image Lab software.

Immunohistochemistry

These experiments were approved by the Human Ethics Committee of Xiangya Hospital, and informed consent was obtained from all patients. Based on polyformalin-fixed and paraffin-embedded tissues obtained from GBM patients, immunohistochemistry analysis was conducted. The tissue sections were first deparaffinized and hydrated for antigen retrieval. They were then incubated with 0.3% hydrogen peroxide for 10 min at room temperature and washed twice with PBS. After blocking with 5% goat serum for 10 min, the sections were washed with PBS and incubated overnight with a primary antibody against the ETV2 antibody (1:100, Abcam, ab181847) at 4°C. A horseradish peroxidase-labeled secondary antibody (1:400, Abcam, United States) was added dropwise to the sections, and incubation was carried out at 37°C for 30 min. After washing with PBS, the sections were developed using a DAB substrate kit (Sangon Biotech, Shanghai, China) and counterstained with hematoxylin (Sangon Biotech).

Statistical Analysis

The R software (version 3.5.1) and GraphPad Prism (version 7.0.0) were used in statistical analyses, and a p < 0.05 was considered significant. The log-rank test was conducted in the Kaplan-Meier survival analysis. The Student’s t-test was used to compare two groups comparing.

Results

The mRNAsi and Clinical Characteristics in Glioma

mRNA expression-based stemness index is a useful indicator that can estimate the number of CSCs by assessing the similarity and heterogeneity between the tumor cells and stem cells. To explore mRNAsi in glioma, the mRNAsi between normal and glioma samples was calculated using TCGA dataset. However, there was no significant difference of mRNAsi between normal and tumor samples (Figure 2A). Next, the correlation between mRNAsi and WHO grades was analyzed. The results showed that different stages of glioma had significantly different mRNAsi values (Figure 2B). Furthermore, the predictive potency of mRNAsi for patient survival prediction was also examined. The Kaplan-Meier analysis showed that higher mRNAsi was significantly associated with a better prognosis of glioma patients (Figure 2C).

FIGURE 2

Figure 2. The mRNA expression-based stemness index (mRNAsi), corrected mRNAsi and clinical characteristics in glioma. (A) The expression level of mRNAsi in normal (five samples) and tumor (606 samples) tissues. (B) Expression level of mRNAsi in gliomas of different WHO grades. (C) Kaplan-Meier survival analysis of the relationship between mRNAsi and overall survival (OS). (D) The expression level of corrected mRNAsi (mRNAsi/tumor purity) in normal (five samples) and tumor (583 samples) tissues (23 samples cannot calculate the tumor purity). (E) The expression level of corrected mRNAsi in glioma of different WHO grades. (F) Kaplan-Meier survival analysis of the relationship between corrected mRNAsi and OS.

Because tumor tissues consist of not only tumor cells but also stromal and immune cells, it reminds us that tumor purity is an important factor interfering with the evaluation of mRNAsi in clinical characteristics. To exclude the potential confounding effect of tumor purity on the analysis, the corrected mRNAsi (mRNAsi/tumor purity) was calculated as previously reported (Pan et al., 2019). We reanalyzed the corrected mRNAsi in normal and glioma samples, but still found no significant difference between them (Figure 2D). However, we found that the corrected mRNAsi values were positively correlated with WHO grades of glioma (Figure 2E). Moreover, we found that patients with higher corrected mRNAsi values had poor prognosis (Figure 2F).

Screening of DEGs and Identification of Key Genes-Related to mRNAsi

First, DEG analysis was performed to compare glioma and normal samples. From this analysis, 3,918 DEGs were screened, of which 1,409 were upregulated, and 2,509 were downregulated (Figure 3A).

FIGURE 3

Figure 3. Screening of critical genes related by mRNAsi. (A) Volcano plot of differentially expressed genes (DEGs); red represents upregulated genes, and green indicates downregulated genes. (B,C) Weighted gene correlation network analysis (WGCNA) of DEGs. Different colors represent different modules. (D) Correlation analysis of the modules and clinical traits with mRNAsi or EREG-mRNAsi. Scatter plot analysis of modules in the red (E) and magenta (F) modules.

To identify the key genes related to mRNAsi, WGCNA analysis was applied based on the selected DEGs. In total, we identified 11 modules (Figures 3B,C), among which the blue and magenta modules exhibited positive correlations with mRNAsi and the red module showed a negative correlation with mRNAsi, with the correlations greater than 0.5 or less than −0.5 (Figure 3D). After calculating the correlations between GS and MM of the three modules, we found that the correlations between GS and MM in the red and magenta modules were highly correlated (more than 0.7; Figures 3E,F). Therefore, the two modules were chosen for further analyses.

Molecular Subtypes of Glioma Based on mRNAsi-Related Genes

In total, the expression profiling of 86 mRNAsi-related genes in 583 samples from TCGA RNA-seq dataset was obtained. The clinical features of these patients are shown in Table 1. Consensus clustering was performed in the 584 samples, and patients could be divided into two significantly different subtypes (S1 and S2 groups; Figures 4A–C). The heatmap of the two subtypes is also been shown in Figure 4D. Compared with patients in the S1 group, glioma patients in the S2 group showed a shorter overall survival (OS; p < 0.0001; Figure 4E). The PPI network analysis showed that most of the mRNAsi-related genes were closely correlated and centered on the stemness-related molecules, such as CD44, CD68, IL6, and CXCR4 (Wang et al., 2019; Ma et al., 2020; Osman et al., 2020).

TABLE 1

Table 1. Clinical characteristics of 583 glioma patients from The Cancer Genome Atlas (TCGA) cohort included in this study.

FIGURE 4

Figure 4. The mRNAsi-related genes could classify glioma into two groups by consensus clustering of TCGA dataset. (A) Cumulative distribution function (CDF) for k = 2 to k = 6. (B) Relative change in area under the CDF curve according to different k values. (C) Consensus clustering matrix of 583 samples from TCGA dataset for k = 2. (D) Heatmap of two clusters defined by the expression of mRNAsi genes. (E) Survival analysis of patients in the S1 group and S2 group in TCGA cohort. (F) Protein-protein interaction (PPI) network of the mRNAsi-related genes.

To validate that mRNAsi-related genes could predict prognostic subtypes, the same method with consensus clustering was applied to the CGGA dataset. Interestingly, the patients also could be divide into two distinct subtypes (S1 and S2 groups; Figures 5A–C). Patients in the S2 group also had significantly worse OS (p < 0.0001; Figure 5D).

FIGURE 5

Figure 5. The mRNAsi-related genes could classify glioma into two groups by consensus clustering of the Chinese Glioma Genome Atlas (CGGA) dataset. (A) CDF for k = 2 to k = 6. (B) Relative change in area under the CDF curve according to different k values. (C) Consensus clustering matrix of 408 samples from TCGA dataset for k = 2. (D) Survival analysis of patients in the S1 group and S2 group in TCGA cohort.

Clinical Characteristics of Two Stemness Subtypes in the Glioma

To identify the clinical characteristics of the two stemness subtypes, we compared the age, IDH status, and WHO grades of the two stemness subtypes in TCGA dataset. There were significantly more elderly patients (age >= 60) in the S2 group than in the S1 group (p < 0.0001; Figure 6A). Moreover, we found that the S1 group had more patients with IDH mutantions (p < 0.0001; Figure 6B). Furthermore, the S1 group had more patients with WHO grade II glioma and fewer patients with WHO IV glioma. However, the S2 group had more patients with WHO grade IV glioma and fewer patients with WHO II glioma (p < 0.0001; Figure 6C). More importantly, these clinical characteristics of the two stemness subtypes could also be validated in the CGGA dataset (Figures 6D–F).

FIGURE 6

Figure 6. Comparison of the clinical characteristics between the two subtypes using TCGA and CGGA datasets. Histograms to showing that the S2 group of TCGA dataset had significantly more elderly patients (A) and more patients with IDH wild-type (B). (C) Pie charts show that the S2 group of TCGA datasets has more patients with WHO grade IV glioma and fewer patients with WHO II glioma than the S1 group (p < 0.001). Histograms to show that the S2 group of CGGA dataset has significantly more elder patients (D) and more patients with IDH wild-type (E). (F) Pie charts show that the S2 group of TCGA datasets has more patients with WHO grade IV glioma and fewer patients with WHO II glioma than the S1 group (p < 0.001). (G) Kaplan-Meier survival curve of S1 and S2 group in different WHO grades (TCGA dataset). (H) Kaplan-Meier survival curve of S1 and S2 group in different WHO grades (TCGA dataset).

Subsequently, we evaluated the prognostic predictive value of the stemness subtypes in different grades. Considering that almost all the WHO IV glioma patients belong to the S2 group, survival analyses were only performed in patients with WHO grade II and III glioma. Stratified survival analyses showed that patients in the S2 group have a better prognosis than patients in the S1 group in both WHO II and III glioma patients in TCGA database (Figure 6G) and CGGA database (Figure 6H).

Immunological Microenvironment in Stemness Subtypes of Glioma

The tumor microenvironment consists of the stromal and immune cells and plays a vital role in the aggressiveness of solid tumors. To measure the level of infiltrating immune cells in the tumor microenvironment, we also used MCP-counter estimates to examine the glioma samples in TCGA database. MCP analysis demonstrated that tumor-associated fibroblasts (CAFs) were significantly higher in the S2 group than in the S1 group (Figures 7A,B). Moreover, there was a similar finding in the CGGA dataset (Figures 7C,D).

FIGURE 7

Figure 7. The immune cell infiltrate distribution of the two stemness subtypes. (A) Heatmap of immune infiltration cells in the two stemness subtypes of TCGA dataset. (B) The proportion of cancer-related fibroblasts in S1 and S2 group of TCGA dataset. ***represents p value of less than 0.001. (C) Heatmap of immune infiltration cells in the two stemness subtypes of CGGA dataset. (D) The proportion of cancer-related fibroblasts in S1 and S2 group of CCGA dataset. ***represents p value of less than 0.001.

Functional Enrichment Analysis Between Two Stemness Subtypes

To explore potential vital molecules and pathways contributing to different subtypes, we performed GO, KEGG, and GSEA analysis between the two stemness subtypes. Figures 8A–C shows the top 20 enriched GO terms in biological processes (BP), cellular components (CC), and molecular functions (MF). GO analysis revealed that immune response and cell adhesion were the main terms involved in BP; plasma membrane and extracellular matrix were significantly enriched in CC; calcium ion binding and channel activity were most enriched in MF (Figures 8A–C). The results of the KEGG pathway analysis showed that cell adhesion and immunological related pathways were mainly enriched (Figure 8D). GSEA showed the significantly enriched hallmark terms, including complement and coagulation cascades, cytokine receptor interaction, intestinal immune network for IgA production, and primary immunodeficiency (Figure 8E).

FIGURE 8

Figure 8. Functional enrichment analysis between two stemness subtypes of TCGA dataset. (A) The top 20 terms of biological processes (BP) in DEGs of two stemness subtypes. (B) The top 20 terms of cellular components (CC). (C) The top 20 terms of molecular functions (MF). (D) The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway in DEGs of two stemness subtypes. (E) The top four significantly enriched pathways of gene set enrichment analysis (GSEA).

Development and Validation of a Stemness Subtype-Related Prognostic Signature

Among the 3,129 candidate DEGs between the two stemness subtypes in TCGA dataset, 3,118 were identified as being independently associated with OS in univariate Cox regression analysis. The top 20 genes were used to perform multivariate Cox analysis. Based on the results of LASSO Cox regression analysis, a stemness subtype-related prognostic signature was developed. The risk score was calculated as follows: [Expression of PTRF × 0.08845 + Expression of ELF4 × 0.01153 + Expression of ELF5 × 0.2005 + Expression of BTN2A2 × 0.06099 + Expression of HMX1 × (−0.04178) + Expression level of FAAH × (−0.13481) + Expression of RGS16 × 0.10462 + Expression of IL4I1 × 0.12804 + Expression of LUZP2 × (−0.08583) + Expression of PLAT × 0.2138 + Expression of ETV2 × 0.15523].

In this prognostic signature, eight genes were negatively associated with OS, and three were positively associated with OS (Figure 9A). Based on this prognostic signature, the risk score for each patient was calculated. According to the median cutoff value of risk scores, all patients were divided into high‐ and low-risk groups in both the training (TCGA dataset) and validation cohorts (CGGA dataset). The distribution of living status and time for each patient in the training and validation cohorts are shown in Figures 9B,C. Patients in the high-risk group had a shorter OS than patients in the low-risk group in the training cohort (p < 0.001; Figure 9D). The time-dependent ROC curve analysis showed that the AUC values of 1, 3, and 5 years were 0.897, 0.892, and 0.826 in the training cohort, respectively (p < 0.001; Figure 9E). Furthermore, glioma samples (both LGGs and GBM) with an IDH1-mutant type have lower risk scores than IDH1wild-type samples, and the risk scores in LGGs with IDH1-mutant and 1p/19q codeletion samples have lower risk scores than IDH1-mutant and 1p/19q non-codeletion samples (Figure 9F). Moreover, we assessed this risk score formula in the CGGA dataset and also found that patients with high-risk scores had poor prognosis in the validation cohort (p < 0.001; Figure 9G). The time-dependent ROC curve analysis showed that the AUC values of this risk score formula at 1, 3, and 5 years were 0.779, 0.858, and 0.822 in the validation cohort, respectively (Figure 9H). The risk scores were also consistent well with the molecular subtypes of glioma (Figure 9I).

FIGURE 9

Figure 9. Construction and validation of the stemness subtype-related risk score model. (A) The 11 genes were selected by least absolute shrinkage and selection operator (LASSO) Cox analysis in TCGA dataset. Risk scores and living status for each patient in the training cohort (B) and the validation cohort (C). (D) Kaplan-Meier curves of the OS of each patient in the training cohort. (E) Time-dependent ROC curve analysis of the risk score model in the training cohort (TCGA dataset). (F) The risk score distributions in different molecular subtypes of glioma (TCGA dataset). (G) Kaplan-Meier curves of the OS of each patient in the validation cohort (CCGA dataset). (H) Time-dependent ROC curve analysis of the risk score model in the validation cohort (CCGA dataset). (I) The risk score distributions in different molecular subtypes of glioma (CGGA dataset). ^***p < 0.001 and ^****p < 0.0001.

Development and Evaluation of the Nomogram

In the univariate Cox analysis, the results showed that the risk score, grade, age, IDH status and 1p/19q status were significantly associated with OS. Then, we performed multivariate Cox regression analyses and the results showed that risk score (HR 1.612; 95% CI 1.258–2.067; p = 1.62E-04), grade (HR 1.488; 95% CI 1.111–1.995; p = 7.76E-03), age (HR 1.035; 95% CI 1.022–1.048; p = 8.85E-08), and 1p/19q status (HR 0.571; 95% CI 0.333–0.980; p = 0.042) were independently related to OS (Figure 10A). Based on the risk score and independent prognostic factors (grade, age, and 1p/19q status) in TCGA dataset, we constructed a nomogram model to predict the prognosis of glioma (Figure 10B). The calibration plot showed that the predicted values of OS at 1-, 3-, and 5-years for glioma patients had a good correlation with the actual values (Figure 10C). Then, the ROC curve analysis of the nomogram also showed a satisfactory evaluation for sensitivity and specificity with a 1-year AUC of 0.909, 3-years AUC of 0.922, and 5-year AUC of 0.874 (Figure 10D).

FIGURE 10

Figure 10. The nomogram based on TCGA dataset for survival prediction in glioma patients. (A) Univariate and multivariate Cox regression analysis of clinical features. (B) Development a nomogram for the quantitative prediction of 1-, 3-, and 5-years survival for LGG patients. (C) The calibration curves for predicting glioma patient 1-, 3-, 5-years survival. (D) The 1-, 3-, and 5-year time-dependent ROC curves of the nomogram.

ETV2 Is Involved in the Migration, Invasion, and EMT Process of Glioma Cells

Among the eight genes that were negatively associated with OS, we tested the expression of ETV2 in clinical tissue samples and glioma cell lines. First, we tested the expression of ETV2 in clinical tissue samples with immunohistochemistry. We found that the expression of ETV2 in clinical patients was correlated with WHO grade (Figure 11A). We also tested the expression of ETV2 in the SHG44 and A172 cell lines with q-PCR, and the expression of ETV2 in glioma cell lines (SHG44 and A172) was higher than that in the normal human cell line (HEB; Figure 11B). Next, we evaluated the effect of ETV2 on glioma cancer cell migration and invasion. Our results revealed that ETV2 knockdown dramatically impaired the cell migration ability of SHG44 and A172 cells relative to the control (Figures 11C,D). In in vitro invasion assays, the invasion potential was obviously suppressed due to the depletion of ETV2 (Figures 11E,F). These findings showed that ETV2 was a significant oncogene associated with the metastatic phenotypes of glioma cells.

FIGURE 11

Figure 11. ETV2 is involved in the migration, invasion, and epithelial-mesenchymal transition (EMT) process of glioma cells. (A) Expression of ETV2 in IHC stained images of glioma tissues compared with adjacent non-tumor tissues was showed. For quantification, three para-cancerous tissues and five glioma sample for tumor groups were counted, mean values with SD are given, p-value is calculated by unpaired t-test, ^**p < 0.01 and ^***p < 0.001. (B) q-PCR showed the relative expression of ETV2 in HEB, SHG44, and A172 cells. (C) Western blot showed the ETV2 siRNA knockdown efficiency in SHG44 and A172 cells. (D) Wound healing assay showed that ETV2 knockdown led to profound impairment in cell migration ability of A172 and SHG44 relative to the control. (E) Invasion assays showed that ETV2 depletion resulted in suppression of A172 and SHG44 cell invasion capability. (F) Vimentin and N-cadherin levels were detected by western blot analysis in A172 and SHG44 cells after ETV2 knockdown with siRNA. For (D,E) mean values with SD form three independent experiments are given. p-value is calculated by unpaired t-test, ^***p < 0.001.

EMT may promote the increased migratory capacity and invasiveness of tumor cells (Zhang et al., 2020b). Therefore, we investigated whether ETV2 mediates EMT in glioma. Considering the reduction in epithelial cells, E-cadherin was poorly expressed in glioma, so we tested the expression of vimentin and N-cadherin to analyze the EMT process in glioma (Liu et al., 2019). The expression of a mesenchymal marker, N-cadherin, decreased after ETV2 knockdown via siRNA. Vimentin is an intermediate filament protein that plays key roles in the integration of cytoskeletal functions and cellular migration. N-cadherin was also decreased after ETV2 knockdown both in both SHG44 and A172 cells. These results indicated that ETV2 is required for the EMT process of glioma.

Discussion

Glioma is the most common and invasive primary brain tumor in adults. Tumor recurrence and treatment resistance are the obstacles to the treatment of glioma. CSCs play essential roles in these processes. Investigating the characteristics of CSCs may facilitate diagnosis, treatment, and prognostic prediction (Chai et al., 2018). In this study, we found that novel molecular subtypes, which based on the stemness index, could effectively predict prognosis in glioma patients. Moreover, the clinical characteristics (age, IDH status, and WHO grades) and tumor microenvironment of the two stemness subtypes are different. This typing could also be validated by the external dataset. Based on the DEGs of the two subtypes, we established a risk score model and a nomogram that could effectively predict the OS of glioma patients. Finally, we selected one gene (ETV2) from the risk score model for experimental validation. To our knowledge, this is the first study to provide a new type of glioma based on the mRNAsi-related genes.

Cancer stemness signatures, which are based on gene expression differences, have been applied to assess the clinical prognosis of some types of tumors (Ben-Porath et al., 2008; Eppert et al., 2011; Pinto et al., 2015). In 2018, Malta et al. (2018) proposed the conception of mRNAsi, which is considered as a more comprehensive index to uncover the characteristics of CSCs. Subsequently, mRNAsi was widely applied to reveal the stemness-related characteristics of different cancers, such as lung cancer (Qin et al., 2020; Zhang et al., 2020a), bladder cancer (Pan et al., 2019), endometrial carcinoma (Liu et al., 2020), medulloblastoma (Lian et al., 2019), and breast cancer (Pei et al., 2020). Moreover, mRNAsi was also used to identify the prognostic biomarkers and therapeutic targets associated with CSC characteristics of glioma (Lvu et al., 2020; Xia et al., 2020). Compared with a previous study, we performed the consensus clustering and identified novel molecular subtypes (S1 and S2 group) of glioma based on mRNAsi-related genes. We also constructed a nomogram model to predict the prognosis of glioma, which has potential for clinical application. More importantly, we performed some in vitro experiments to confirm the function of ETV2 in glioma cells.

Gliomas are the most frequent intrinsic tumors of the central nervous system. In the revised fourth edition of the WHO classification of CNS tumors published in 2016, the status of IDH and 1p/19q codeletion was applied in classification (Zhang et al., 2013). In our study, patients in subtype I (S1 group) were younger, more likely to have IDH-mutant status, lower WHO grades and poor prognosis than those in the subtype II (S2 group). Based on our risk score model, we found that glioma samples (both LGGs and GBM) with an IDH1-mutant type have lower risk scores than IDH1wild-type samples, and the risk scores in LGGs with IDH1-mutant and 1p/19q codeletion samples have lower risk scores compared with than IDH1-mutant and 1p/19q non-codeletion samples. The immunological microenvironments between the two subtypes were also different. Gliomas in subtype II (S2 group) are more likely to have a higher proportion of tumor-related immune and stromal cells, especially CAFs. CAFs, one of the main cellular components of the tumor microenvironment, play an important role in promoting cancer cell invasion and dissemination (Hurtado et al., 2020). These results indicated that this new typing could provide novel mechanistic and clinical insights for the diagnosis, treatment, and prognostic prediction of glioma patients.

Functional enrichment analysis of the DEGs between two stemness subtypes showed that immune and ECM-related GO terms and pathways were mainly enriched. This implied that the immunological microenvironment and ECM components might have a close relationship with the stemness characteristics of glioma. Previous studies have shown that glioma stem cells interact with immune cells and simulate the early microenvironment during tumorigenesis (Zhai et al., 2020). CSCs can protect cancer cells from immune attack by producing immune inhibitory factors to communicate with tumor microenvironment components (Khosravi et al., 2020). Moreover, CSCs contribute to glioma invasiveness, which is closely correlated with the extracellular matrix (Ortensi et al., 2013). These findings were in accordance with our findings.

In the present study, we also developed a stemness subtype-related prognostic signature. There were 11 essential prognostic genes in this signature, which have not been reported in previous glioma stemness index related publications (Lvu et al., 2020; Xia et al., 2020). Our results are a good supplement to the existing research about the prognosis prediction of patients with glioma. Among the three prognostic genes that were positively associated with OS, LUZP2 is a protein limitedly expressed in the brain and spinal cord. A recent study showed that low LUZP2 expression independently predicted poor OS in LGG (Li et al., 2020). Fatty-acid amide hydrolase (FAAH), an intracellular serine hydrolase, plays an important role in the inhibition of stem cell migration (Wollank et al., 2015). Among the eight genes that were negatively associated with OS, ETV2, ELF5, IL4I1, and BTN2A2 were novel prognostic biomarkers for glioma, which have not yet been reported yet. ETV2 is a critical factor for vascular development and regeneration, which may contribute to tumorigenesis (Baltrunaite et al., 2017; Choi, 2018). In GBM, ETV2 is sufficient and necessary for the trans-differentiation of GBM stem cells to an endothelial lineage (Humm and Sylvia, 1965). In our study, we found that ETV2 is negatively associated with OS. The expression of ETV2 is closely associated with WHO grade. More importantly, in vitro experiment revealed that ETV2 is involved in the migration, invasion, and EMT process of glioma cells. ELF5 is an epithelial-specific member of the E26 transforming sequence (ETS) transcription factor family, which plays critical roles in malignancy, particularly in basal-like and endocrine-resistant forms of breast cancer (Piggin et al., 2016; Singh et al., 2020). IL4I1 and BTN2A2 are both involved in the regulation of the immunologic microenvironment in different tumors (Smith et al., 2010; Sarter et al., 2016; Molinier-Frenkel et al., 2019; Sadik et al., 2020). However, what is their roles in tumorigenesis and progression of glioma are still need further investigation.

There were some limitations in our study. First, there were only very limited normal samples (only five normal samples) included in our study, which might lead to there being no difference in mRNAsi between normal and tumor samples. Second, DEG analysis might neglect some potential mRNAs that were closely related to the mRNAsi. Third, this study is a retrospective study, and the stemness-related typing of glioma should be further confirmed by prospective studies. Finally, the underlying mechanisms of the selected genes in the model affecting the prognosis of glioma should be elucidated by more experiments in vivo.

In this study, we identified two distinct stemness-related molecular subtypes of glioma, which could provide new insights for the development of precision diagnosis and prognostic prediction for glioma patients. Moreover, we developed a stemness subtype-related prognostic signature that could effectively predict the prognosis of glioma patients.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author Contributions

XJ and WY conceived and designed the study. JT and WY wrote the manuscript. YC, ZX, QZ, and CZ analyzed the results. GT, HL, SW, YG, and ZJ performed the image visualization. JT, HZ, MZ, and CR performed the validation experiment. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (no. 81472355) and the Hunan Provincial Science and Technology Department (no. 2014FJ6006).

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.

Footnotes

1. https://xena.ucsc.edu/public

2. http://www.cgga.org.cn/index.jsp

3. https://string-db.org/

4. http://software.broadinstitute.org/gsea/msigdb/

References

Baltrunaite, K., Craig, M. P., Palencia Desai, S., Chaturvedi, P., Pandey, R. N., Hegde, R. S., et al. (2017). ETS transcription factors Etv2 and Fli1b are required for tumor angiogenesis. Angiogenesis 20, 307–323. doi: 10.1007/s10456-017-9539-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17:218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Porath, I., Thomson, M. W., Carey, V. J., Ge, R., Bell, G. W., Regev, A., et al. (2008). An embryonic stem cell-like gene expression signature in poorly differentiated aggressive human tumors. Nat. Genet. 40, 499–507. doi: 10.1038/ng.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Bush, N. A., Chang, S. M., and Berger, M. S. (2017). Current and future strategies for treatment of glioma. Neurosurg. Rev. 40, 1–14. doi: 10.1007/s10143-016-0709-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chai, R., Zhang, K., Wang, K., Li, G., Huang, R., Zhao, Z., et al. (2018). A novel gene signature based on five glioblastoma stem-like cell relevant genes predicts the survival of primary glioblastoma. J. Cancer Res. Clin. Oncol. 144, 439–447. doi: 10.1007/s00432-017-2572-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, K. (2018). ETS transcription factor ETV2/ER71/Etsrp in haematopoietic regeneration. Curr. Opin. Hematol. 25, 253–258. doi: 10.1097/MOH.0000000000000430

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, M. F. (2019). Clinical and therapeutic implications of cancer stem cells. N. Engl. J. Med. 380, 2237–2245. doi: 10.1056/NEJMra1804280

PubMed Abstract | CrossRef Full Text | Google Scholar

Eppert, K., Takenaka, K., Lechman, E. R., Waldron, L., Nilsson, B., van Galen, P., et al. (2011). Stem cell gene expression programs influence clinical outcome in human leukemia. Nat. Med. 17, 1086–1093. doi: 10.1038/nm.2415

PubMed Abstract | CrossRef Full Text | Google Scholar

Goeman, J. J. (2010). L1 penalized estimation in the cox proportional hazards model. Biom. J. 52, 70–84. doi: 10.1002/bimj.200900028

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, W. D., Sherman, B. T., and Lempicki, R. A. (2009a). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, W. D., Sherman, B. T., and Lempicki, R. A. (2009b). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37, 1–13. doi: 10.1093/nar/gkn923

PubMed Abstract | CrossRef Full Text | Google Scholar

Humm, D. G., and Sylvia, A. L. (1965). Soluble proteins of a melanoma and normal skin from the swordtail, platyfish, and their hybrids. Science 150, 635–636. doi: 10.1126/science.150.3696.635

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurtado, P., Martínez-Pena, I., and Piñeiro, R. (2020). Dangerous liaisons: circulating tumor cells (CTCs) and cancer-associated fibroblasts (CAFs). Cancer 12:2861. doi: 10.3390/cancers12102861

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, W., Jiang, X., Liu, W., Wang, L., Zhu, B., Zhu, H., et al. (2018). Effects of three-dimensional collagen scaffolds on the expression profiles and biological functions of glioma cells. Int. J. Oncol. 52, 1787–1800. doi: 10.3892/ijo.2018.4330

PubMed Abstract | CrossRef Full Text | Google Scholar

Khosravi, N., Mokhtarzadeh, A., Baghbanzadeh, A., Hajiasgharzadeh, K., Shahgoli, V. K., Hemmat, N., et al. (2020). Immune checkpoints in tumor microenvironment and their relevance to the development of cancer stem cells. Life Sci. 256:118005. doi: 10.1016/j.lfs.2020.118005

PubMed Abstract | CrossRef Full Text | Google Scholar

Klughammer, J., Kiesel, B., Roetzer, T., Fortelny, N., Nemc, A., Nenning, K. H., et al. (2018). The DNA methylation landscape of glioblastoma disease progression shows extensive heterogeneity in time and space. Nat. Med. 24, 1611–1624. doi: 10.1038/s41591-018-0156-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Deng, G., Qi, Y., Zhang, H., Jiang, H., Geng, R., et al. (2020). Downregulation of LUZP2 Is correlated with poor prognosis of low-grade glioma. Biomed. Res. Int. 2020:9716720. doi: 10.1155/2020/9716720

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, H., Han, Y. P., Zhang, Y. C., Zhao, Y., Yan, S., Li, Q. F., et al. (2019). Integrative analysis of gene expression and DNA methylation through one-class logistic regression machine learning identifies stemness features in medulloblastoma. Mol. Oncol. 13, 2227–2245. doi: 10.1002/1878-0261.12557

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, N., Wang, Z., Liu, D., and Xie, P. (2019). HOXC13-AS-miR-122-5p-SATB1-C-Myc feedback loop promotes migration, invasion and EMT process in glioma. OncoTargets Ther. 12, 7165–7173. doi: 10.2147/OTT.S220027

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Wu, Z., Sun, R., Nie, S., Meng, H., Zhong, Y., et al. (2020). Using mRNAsi to identify prognostic-related genes in endometrial carcinoma based on WGCNA. Life Sci. 258:118231. doi: 10.1016/j.lfs.2020.118231

PubMed Abstract | CrossRef Full Text | Google Scholar

Lvu, W., Fei, X., Chen, C., and Zhang, B. (2020). In silico identification of the prognostic biomarkers and therapeutic targets associated with cancer stem cell characteristics of glioma. Biosci. Rep. 40:BSR20201037. doi: 10.1042/BSR20201037

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, X., Chen, J., Liu, J., Xu, B., Liang, X., Yang, X., et al. (2020). IL-8/CXCR2 mediates tropism of human bone marrow-derived mesenchymal stem cells toward CD133(+) /CD44(+) colon cancer stem cells. J. Cell. Physiol. 236, 3114–3128. doi: 10.1002/jcp.30080

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell 173, 338.e15–354.e15. doi: 10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Molinier-Frenkel, V., Prévost-Blondel, A., and Castellano, F. (2019). The IL4I1 enzyme: a new player in the immunosuppressive tumor microenvironment. Cells 8:757. doi: 10.3390/cells8070757

PubMed Abstract | CrossRef Full Text | Google Scholar

Nabors, L. B., Portnow, J., Ammirati, M., Baehring, J., Brem, H., Butowski, N., et al. (2017). NCCN guidelines insights: central nervous system cancers, version 1.2017. J. Natl. Compr. Cancer Netw. 15, 1331–1345. doi: 10.6004/jnccn.2017.0166

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortensi, B., Setti, M., Osti, D., and Pelicci, G. (2013). Cancer stem cell contribution to glioblastoma invasiveness. Stem Cell Res. Ther. 4:18. doi: 10.1186/scrt166

PubMed Abstract | CrossRef Full Text | Google Scholar

Osman, A., Oze, M., Afify, S. M., Hassan, G., El-Ghlban, S., Nawara, H. M., et al. (2020). Tumor-associated macrophages derived from cancer stem cells. Acta Histochem. 122:151628. doi: 10.1016/j.acthis.2020.151628

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostrom, Q. T., Bauchet, L., Davis, F. G., Deltour, I., Fisher, J. L., Langer, C. E., et al. (2014). The epidemiology of glioma in adults: a "state of the science" review. Neuro Oncol. 16, 896–913. doi: 10.1093/neuonc/nou087

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, S., Zhan, Y., Chen, X., Wu, B., and Liu, B. (2019). Identification of biomarkers for controlling cancer stem cell characteristics in bladder cancer by network analysis of transcriptome data stemness indices. Front. Oncol. 9:613. doi: 10.3389/fonc.2019.00613

PubMed Abstract | CrossRef Full Text | Google Scholar

Pei, J., Wang, Y., and Li, Y. (2020). Identification of key genes controlling breast cancer stem cell characteristics via stemness indices analysis. J. Transl. Med. 18:74. doi: 10.1186/s12967-020-02260-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Piggin, C. L., Roden, D. L., Gallego-Ortega, D., Lee, H. J., Oakes, S. R., and Ormandy, C. J. (2016). ELF5 isoform expression is tissue-specific and significantly altered in cancer. Breast Cancer Res. 18:4. doi: 10.1186/s13058-015-0666-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, J. P., Kalathur, R. K., Oliveira, D. V., Barata, T., Machado, R. S., Machado, S., et al. (2015). StemChecker: a web-based tool to discover and explore stemness signatures in gene sets. Nucleic Acids Res. 43, W72–W77. doi: 10.1093/nar/gkv529

PubMed Abstract | CrossRef Full Text | Google Scholar

Prieto-Vila, M., Takahashi, R. U., Usuba, W., Kohama, I., and Ochiya, T. (2017). Drug resistance driven by cancer stem cells and their niche. Int. J. Mol. Sci. 18:2574. doi: 10.3390/ijms18122574

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, S., Long, X., Zhao, Q., and Zhao, W. (2020). Co-expression network analysis identified genes associated with cancer stem cell characteristics in lung squamous cell carcinoma. Cancer Invest. 38, 13–22. doi: 10.1080/07357907.2019.1697281

PubMed Abstract | CrossRef Full Text | Google Scholar

Sadik, A., Somarribas Patterson, L. F., Öztürk, S., Mohapatra, S. R., Panitz, V., Secker, P. F., et al. (2020). IL4I1 is a metabolic immune checkpoint that activates the AHR and promotes tumor progression. Cell 182, 1252.e34–1270.e34. doi: 10.1016/j.cell.2020.07.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarter, K., Leimgruber, E., Gobet, F., Agrawal, V., Dunand-Sauthier, I., Barras, E., et al. (2016). Btn2a2, a T cell immunomodulatory molecule coregulated with MHC class II genes. J. Exp. Med. 213, 177–187. doi: 10.1084/jem.20150435

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz, A., Meyer, F., Dubrovska, A., and Borgmann, K. (2019). Cancer stem cells and radioresistance: DNA repair and beyond. Cancer 11:862. doi: 10.3390/cancers11060862

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, S., Kumar, S., Srivastava, R. K., Nandi, A., Thacker, G., Murali, H., et al. (2020). Loss of ELF5-FBXW7 stabilizes IFNGR1 to promote the growth and metastasis of triple-negative breast cancer through interferon-γ signalling. Nat. Cell Biol. 22, 591–602. doi: 10.1038/s41556-020-0495-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, I. A., Knezevic, B. R., Ammann, J. U., Rhodes, D. A., Aw, D., Palmer, D. B., et al. (2010). BTN1A1, the mammary gland butyrophilin, and BTN2A2 are both inhibitors of T cell activation. J. Immunol. 184, 3514–3525. doi: 10.4049/jimmunol.0900416

PubMed Abstract | CrossRef Full Text | Google Scholar

Vlashi, E., and Pajonk, F. (2015). Cancer stem cells, cancer cell plasticity and radiation therapy. Semin. Cancer Biol. 31, 28–35. doi: 10.1016/j.semcancer.2014.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Li, X., Wu, H., Wang, H., Yao, L., Deng, Z., et al. (2019). EMP1 regulates cell proliferation, migration, and stemness in gliomas through PI3K-AKT signaling and CD44. J. Cell. Biochem. 120, 17142–17150. doi: 10.1002/jcb.28974

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wollank, Y., Ramer, R., Ivanov, I., Salamon, A., Peters, K., and Hinz, B. (2015). Inhibition of FAAH confers increased stem cell migration via PPARα. J. Lipid Res. 56, 1947–1960. doi: 10.1194/jlr.M061473

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, P., Li, Q., Wu, G., and Huang, Y. (2020). Identification of glioma cancer stem cell characteristics based on weighted gene prognosis module co-expression network analysis of transcriptome data stemness indices. J. Mol. Neurosci. 70, 1512–1520. doi: 10.1007/s12031-020-01590-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, W., Jiang, X., Tan, J., Xin, Z., Zhou, Q., Zhan, C., et al. (2020). Development and validation of a tumor mutation burden-related immune prognostic model for lower-grade glioma. Front. Oncol. 10:1409. doi: 10.3389/fonc.2020.01409

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhai, Y., Li, G., Li, R., Chang, Y., Feng, Y., Wang, D., et al. (2020). Single-cell RNA-sequencing shift in the interaction pattern between glioma stem cells and immune cells during tumorigenesis. Front. Immunol. 11:581209. doi: 10.3389/fimmu.2020.581209

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Bao, Z., Zhang, W., and Jiang, T. (2013). Progress on molecular biomarkers and classification of malignant gliomas. Front. Med. 7, 150–156. doi: 10.1007/s11684-013-0267-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Tseng, J. T., Lien, I. C., Li, F., Wu, W., and Li, H. (2020a). mRNAsi index: machine learning in mining lung adenocarcinoma stem cell biomarkers. Genes 11:257. doi: 10.3390/genes11121395

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Yang, L., Guan, G., Cheng, P., Cheng, W., and Wu, A. (2020b). LOXL2 upregulation in gliomas drives tumorigenicity by activating autophagy to promote TMZ resistance and trigger EMT. Front. Oncol. 10:569584. doi: 10.3389/fonc.2020.569584

PubMed Abstract | CrossRef Full Text | Google Scholar