key: cord-0012018-b9zlig0r
authors: Mukherjee, Shradha
title: Quiescent stem cell marker genes in glioma gene networks are sufficient to distinguish between normal and glioblastoma (GBM) samples
date: 2020-07-02
journal: Sci Rep
DOI: 10.1038/s41598-020-67753-5
sha: f4b1a85ef5cef5d8512e331283984d6421d3d942
doc_id: 12018
cord_uid: b9zlig0r

Grade 4 glioma or GBM has poor prognosis and is the most aggressive grade of glioma. Accurate diagnosis and classification of tumor grade is a critical determinant for development of treatment pathway. Extensive genomic sequencing of gliomas, different cell types, brain tissue regions and advances in bioinformatics algorithms, have presented an opportunity to identify molecular markers that can complement existing histology and imaging methods used to diagnose and classify gliomas. ‘Cancer stem cell theory’ purports that a minor population of stem cells among the heterogeneous population of different cell types in the tumor, drive tumor growth and resistance to therapies. However, characterization of stem cell states in GBM and ability of stem cell state signature genes to serve as diagnostic or prognostic molecular markers are unknown. In this work, two different network construction algorithms, Weighted correlation network analysis (WGCNA) and Multiscale Clustering of Geometric Network (MEGENA), were applied on publicly available glioma, control brain and stem cell gene expression RNA-seq datasets, to identify gene network regulatory modules associated with GBM. Both gene network algorithms identified consensus or equivalent modules, HuAgeGBsplit_18 (WGCNA) and c1_HuAgeGBsplit_32/193 (MEGENA), significantly associated with GBM. Characterization of HuAgeGBsplit_18 (WGCNA) and c1_HuAgeGBsplit_32/193 (MEGENA) modules showed significant enrichment of rodent quiescent stem cell marker genes (GSE70696_QNPbyTAP). A logistic regression model built with eight of these quiescent stem cell marker genes (GSE70696_QNPbyTAP) was sufficient to distinguish between control and GBM samples. This study demonstrates that GBM associated gene regulatory modules are characterized by diagnostic quiescent stem cell marker genes, which may potentially be used clinically as diagnostic markers and therapeutic targets in GBM.

SVA + LM adjustment the resultant normalized log2TPM + 1 gene expression was used to build a glioma gene network with subnetworks or modules using WGCNA_1.66 and MEGENA_1.3.7 R packages 23, 24, 35, 46 .

In WGCNA, default parameters and a minimum module size of 100 was used to calculate a topology overlap matrix (TOM) based on gene expression correlations. Hierarchical clustering was then used to build a glioma gene network consisting of interconnected subnetworks or modules 22, 23 . In MEGENA, default parameters and a minimum module size of 100 was used to calculate a planar filtered network (PFN) from gene expression correlations. Multiscale clustering method was applied to build a glioma gene network consisting of interconnected subnetworks or modules 24 . WGCNA computes scale-free or single scale networks, while MEGENA computes multi-scale networks to include different possible variations of gene-gene interactions. Therefore, in MEGENA a given gene or node can be part of multiple modules representing different possible interactions, while in WGCNA a given gene or node is assigned to only a single module. To determine module trait correlations, module eigengenes were computed with moduleEigengenes R function and correlations were calculated 22 . To compare WGCNA and MEGENA modules, previously published module preservation analysis and hypergeometric enrichment tests were used 35, 47 . Briefly, hypergeometric test was implemented with userListEnichment R function widely used to compare WGCNA gene network modules with each other and with user supplied gene lists 22, 23, 46 . Both module preservation and userListEnrichment R functions are from WGCNA 23, 47 . After module preservation analysis between WGCNA and MEGENA modules, significant overlap of genes between WGCNA and MEGENA modules was done with userListEnrichment R function for all WGCNA and MEGENA modules i.e. each WGCNA module was compared with each MEGENA module. WGCNA and MEGENA modules that significantly overlapped and were significantly associated with GBM were retained.

Stem cell differential gene expression analysis with limma, edgeR and simple comparison of means. Differentially expressed genes (DEGs) between proliferative and quiescent stem cell states were identified using R packages limma_3.38.3, edgeR_3.24.3 and simple comparison expression means 18, 45, 48 . To calculate genes enriched in proliferative stem cells, gene expression of samples annotated to proliferative stem cells were compared with gene expression of samples annotated to quiescent stem cells for each of the five datasets (GSE68270, GSE70696, GSE99777, GSE93991 and GSE114574). Similarly, to calculate genes enriched in quiescent stem cells, gene expression of samples annotated to quiescent stem cells were compared with gene expression of samples annotated to proliferative stem cells for each of the five datasets (GSE68270, GSE70696, GSE99777, GSE93991 and GSE114574). In simple comparison of means method, mean expression of genes were simply compared between proliferative and quiescent stem cell states to determine DEGs. In limma and edgeR model design included variables stem cell state, study or batch, gender, age and tissue, and stem cell states were contrasted to determine DEGs, while all other variables were held constant in the model. In limma and edgeR methods Benjamini and Hochberf (BH) corrections for multiple testing is included as a large number of genes were included in analysis 49 . DEGs with BH corrected adjP-values < 0.05 and fold change > 1. 25 were considered significant DEGs. To visualize gene expression values of DEGs volcano plots, barplots and density plots were made using ggplot2_3.1.1 R package 40 . Consensus DEGs were obtained by overlapping DEG lists produced by limma, edgeR and simple comparison of means with a significance of overlap p-value < 0.05 as calculated with GeneOverlap_1.18.0 and visualized by VennDiagram_1.6.20 R packages 50, 51 . Consensus DEGs that belonged to atleast two of the three DEG lists produced by limma, edgeR and simple comparison of means were designated significantly enriched genes or DEGs in proliferative and quiescent stem cell states-simply referred to as (1) proliferative stem cell marker genes and (2) quiescent stem cell marker genes.

Following is a detailed description of all DEG analysis contrasts, sample size for each dataset and abbreviations used to represent proliferative and quiescent stem cell marker genes: (A) Adult proliferative neural progenitor cells (PNPCs) vs adult quiescent neural stem cells (QNSCs) DEG analysis to identify genes enriched in PNPCs relative to QNSCs in mouse dataset with series number GSE68270 and sample size of 4 each, abbreviated as GSE68270_PNPCbyQNSC 32 (B) Adult quiescent neural stem cells (QNSCs) vs adult proliferative neural progenitor cells (PNPCs) DEG analysis to identify genes enriched in QNSCs relative to PNPCs in mouse dataset with series number GSE68270 and sample size of 4 each, abbreviated as GSE68270_QNSCbyPNPC 32 (C) Adult hippocampal stem cells in proliferative condition or transient amplifying progenitor cells (TAPs) vs adult hippocampal stem cells in quiescent condition or quiescent progenitor cells (QNPs) DEG analysis to identify genes enriched in TAPs relative to QNPs in rat dataset with series number GSE70696 and sample size of 2 each, abbreviated as GSE70696_TAPbyQNP 33 (D) Adult hippocampal stem cells in quiescent condition or quiescent progenitor cells (QNPs) vs adult hippocampal stem cells in proliferative condition or transient amplifying progenitor cells (TAPs) DEG analysis to identify genes enriched in QNPs relative to TAPs in rat dataset with series number GSE70696 and sample size of 2 each, abbreviated as GSE70696_QNPbyTAP 33 (E) Adult proliferative sub-ventricular zone stem cells (PSVZSCs) vs adult quiescent sub-ventricular zone stem cells (QSVZSCs) DEG analysis to identify genes enriched in PSVZSCs relative to QSVZSCs in mouse microarray dataset with series number GSE99777 and sample size of 3 each, abbreviated as GSE99777_PSVZSCbyQSVZSC 34 (F) Adult quiescent sub-ventricular zone stem cells (QSVZSCs) vs adult proliferative sub-ventricular zone stem cells (PSVZSCs) DEG analysis to identify genes enriched in QSVZSCs relative to PSVZSCs in mouse microarray dataset with series number GSE99777 and sample size of 3 each, abbreviated as GSE99777_QSVZSCbyPSVZSC 34 (G) GBM cells cultured in proliferative condition or proliferative GBM cells (PGBCs) vs GBM cells cultured in quiescent condition or quiescent GBM cells (QGBCs) DEG analysis to identify genes enriched in PGBCs relative to QGBCs in human dataset with series number GSE93991 and sample size of 9 and 6, respectively, abbreviated as GSE93991_PGBCbyQGBC 30 (H) GBM cells cultured in quiescent condition or quiescent GBM cells (QGBCs) vs GBM cells cultured in proliferative condition or proliferative GBM cells (PGBCs) DEG analysis to identify genes enriched in QGBCs relative to PGBCs in human dataset with series number GSE93991 and sample size of 6 and Scientific RepoRtS | (2020) 10:10937 | https://doi.org/10.1038/s41598-020-67753-5 www.nature.com/scientificreports/ 9, respectively, abbreviated as GSE93991_QGBCbyPGBC 30 (I) GBM organoids cultured in proliferative condition or proliferative GBM organoids (PGBOs) vs GBM organoids cultured in quiescent condition or quiescent GBM cells (QGBOs) DEG analysis to identify genes enriched in PGBOs relative to QGBOs in human dataset with series number GSE114574 and sample size 6, abbreviated as GSE114574_PGBObyQGBO 31 and (J) GBM organoids cultured in quiescent condition or quiescent GBM organoids (QGBOs) vs GBM organoids cultured in proliferative condition or proliferative GBM organoids (PGBOs) DEG analysis to identify genes enriched in QGBOs relative to PGBOs in human dataset with series number GSE93991 and sample size of 6, abbreviated as GSE114574_QGBObyPGBO 31 .

Enrichment of proliferative and quiescent stem cell marker genes in glioma modules. Enrichment of proliferative and quiescent stem cell marker genes identified by differential gene expression analysis above, in WGCNA and MEGENA modules was determined using useListEnrichment R function 22, 23, 46 . Additionally, a supplementary table containing a set of 336 genes potentially involved in transition of GBM from stem-like state to differentiation identified by SWIM tool were downloaded directly from the published paper 52, 53 . Enrichment of SWIM GBM list in WGCNA and MEGENA modules was also determined using useListEnrichment R function 22, 23, 46 

To determine proliferative and quiescent stem cell marker genes common between all five stem cell datasets (GSE68270, GSE70696, GSE99777, GSE93991 and GSE114574), the following consensus significantly enriched genes or DEGs in proliferative and quiescent states were overlapped, respectively, using online tool https ://www.molbi otool s.com/listc ompar e.html : (A) For proliferation comparison of (GSE68270_PNPCbyQNSC + GSE70696_ TAPbyQNP + GSE99777_PSVZSCbyQSVZSC + GSE93991_PGBCbyQGBC + GSE114574_PGBObyQGBO) and (B) For quiescence comparison of (GSE68270_QNSCbyPNPC + GSE70696_QNPbyTAP + GSE99777_ QSVZSCbyPSVZSC + GSE93991_QGBCbyPGBC + GSE114574_QGBObyPGBO).

To determine proliferative and quiescent stem cell marker genes common to both normal stem cells and GBM in culture, following consensus significantly enriched genes or DEG lists in proliferative and quiescent states were overlapped, respectively, using online tool https ://www.molbi otool s.com/listc ompar e.html (A) comparison of normal stem cells (GSE68270_PNPCbyQNSC + GSE70696_TAPbyQNP + GSE99777_ PSVZSCbyQSVZSC) and GBM cell cultures (GSE93991_PGBCbyQGBC + GSE114574_PGBObyQGBO) DEG 

SVA + LM approach reduces batch effects in RNA-seq datasets. To identify GBM specific transcriptome features, a meta-analysis was performed with glioma and control brain human RNA-seq samples. SVA + LM normalization reduced variability due to batch effects as indicated by greater overlap of expression data in density plots after SVA + LM normalization (Fig. 1A) . Box-Whiskers plots showed a slightly skewed mean expression in GSE67333 dataset that was fixed by SVA + LM normalization (Fig. 1B) . PCA plots showed that SVA + LM normalization reduces dispersion of samples from same study (Fig. 1C) . Correlation plots to evaluate correlation of gene expression values among different studies, showed an increase in positive correlation after SVA + LM normalization (Fig. 1D) . Thus, SVA + LM normalization may be used to achieve reduction in batch effects in global gene expression when RNA-seq studies are combined.

Network analysis reveals GBM associated modules within glioma network. Glioma gene coexpression networks were constructed with WGCNA and MEGENA to uncover underlying molecular mechanisms. WGCNA identified 39 modules with largest module HuAgeGBsplit_01 comprising of 7,969 genes and smallest module HuAgeGBsplit_38 comprising of 113 genes (Table 1) . MEGENA identified 235 modules with largest module c1_HuAgeGBsplit_13 comprising of 1,483 genes and smallest modules c1_HuAgeGB-split_1081/181/1895 comprising of 100 genes each ( Table 2) .

To identify GBM associated modules, a module eigengene was used to represent overall expression pattern of each module produced by WGCNA and MEGENA. Spearman correlations were calculated for GBM and other clinical traits, such as batch, age, and gender. Eight modules in WGCNA, with module size ranging from 673 genes in HuAgeGBsplit_08 to 121 genes in Hu_AgeGBsplit_36, were found to be significantly associated with GBM ( Fig. 2A,B) . All WGCNA modules, including GBM associated modules, were preserved with MEGENA modules (Fig. 2C) . Comparison of GBM specific WGCNA modules with MEGENA modules showed a significant overlap of all WGCNA GBM modules with MEGENA modules (Fig. 2D ). All 20 MEGENA modules that significantly overlapped with GBM WGCNA modules, with module size ranging from 101 genes in c1_HuAgeGBsplit_605 to 708 genes in c1_HuAgeGBsplit_24, were also significantly associated with GBM ( Fig. 2E,F) . Thus, WGCNA and MEGENA complemented each other and helped identify GBM specific modules in largely preserved glioma WGCNA and MEGENA gene networks.

Differential gene expression analysis reveals proliferative and quiescent stem cell marker genes. To identify genes specific to quiescent and proliferative states of stem cells, differential gene expression analysis was performed on different stem cell datasets as described under methods. In mouse adult hippocampal stem cell dataset, 5,306 (3,959 human gene symbols) and 4,072 (3,188 human gene symbols) DEGs with a fold change of 1.25 (p-value < 0.05) were identified in proliferation (GSE68270_PNPCbyQNSC) and quiescence (GSE68270_PNPCbyQNSC), respectively (Fig. 3A,B) . In rat adult hippocampal stem cell dataset, 5,122 (4,362 human gene symbols) and 3,290 (2,792 human gene symbols) DEGs with a fold change of 1.25 (p-value < 0.05) were identified in proliferation (GSE70696_TAPbyQNP) and quiescence (GSE70696_QNPbyTAP), respectively ( Fig. 3C,D) . In mouse adult subventricular zone stem cell dataset, 5,216 (4,235 human gene symbols) and 5,733 (4,658 human gene symbols) DEGs with a fold change of 1.25 (p-value < 0.05) were identified in proliferation (GSE99777_PSVZSCbyQSVZSC) and quiescence (GSE99777_QSVZSCbyPSVZSC), respectively (Fig. 3E,F) . In human GBM cell culture dataset, 7,857 and 3,363 (human gene symbols) DEGs with a fold change of 1.25 (p-value < 0.05) were identified in proliferation (GSE93991_PGBCbyQGBC) and quiescence (GSE93991_QGB-CbyPGBC), respectively (Fig. 3G,H) . In human GBM organoid culture dataset, 1928 and 3,315 (human gene symbols) DEGs with a fold change of 1.25 (p-value < 0.05) were identified in proliferation (GSE114574_PGBObyQGBO) and quiescence (GSE114574_QGBObyPGBO), respectively (Fig. 3I ,J). Stem cell marker genes common between all normal stem cells and GBM culture datasets (GSE68270, GSE70696, GSE99777, GSE93991 and GSE114574), consisted of 7 proliferation genes or DEGs (ACYP1, AKAP12, LRP11, MYSM1, SLC20A1, TERT, TSPAN13) and 5 quiescence genes or DEGs (ETHE1, FZD9, NINJ1, P2RX4, PTP4A3) ( Tables 3, 4 ). Overlapping stem cell marker genes obtained from normal stem cell datasets (GSE68270, GSE70696 and GSE99777), with GBM culture datasets (GSE93991 and GSE114574), showed an overlap of 4,176 proliferation genes (45.87% of GBM proliferation DEGs) and only 1598 quiescence genes (25.97% of GBM quiescence DEGs (Supplementary Tables S1A,B ).

showed that HuAgeGBsplit_18 WGCNA GBM module and its equivalent c1_HuAgeGBsplit_193/32 MEGENA GBM modules, were the only GBM modules significantly enriched with proliferative and quiescent stem cell marker genes (Fig. 4A,B) . Comparison of genes in c1_HuAgeGBsplit_193 and c1_HuAgeGBsplit_32 MEGENA modules, revealed that all genes in c1_HuAgeGBsplit_193 were also present in c1_HuAgeGBsplit_32 ( 

Grade_3 Grade_4 Grade_4 CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CON CONCON CON CON CON CON CON CON CON CON CON CON 

A different color is used for each study Logistic regression model built with quiescent stem cell marker genes in GBM modules. A logistic regression model was built with select quiescent stem cell marker genes (GSE70696_QNPbyTAP) to diagnostic between control and GBM samples. From a total of 110 genes from GSE70696_QNPbyTAP enriched in GBM WGCNA module HuAgeGBsplit_18, genes that were atleast 40-fold upregulated in QNP relative to TAP were selected (CD151, CEND1, DCHS1, SMPD1, TPP1, GATD1, RNH1 and SMCR8) for logistic regression. Effects plots showed that probability of GBM relative to control increases with increased expression of CEND1, DCHS1, TPP1, GATD1, RNH1 and SMCR8 (Fig. 5B ,C,E-H) and decreases with increased expression of CD151 and SMPD1 (Fig. 5A,D) . Logistic regression model with these 8 genes without gene-gene interaction term was not significant (Chi-square p-value = 0.9847), while with gene-gene interactions the model was a significant (Chi-square p-value = 0.00799) predictor for GBM (Fig. 5I ,J).

Hosmer-Lemeshow GOF test on without gene-gene interaction logistic regression model showed a large difference between observed and expected probabilities, and there was significant evidence of poor fit (p-value = 0.006, less than 0.05) (Fig. 5I , Table 5 ). Therefore, in logistic regression model without gene-gene interaction, Ho is rejected and Ha is accepted, and the model is rejected for being a poor fit for the data. Hosmer-Lemeshow GOF test on with gene-gene interaction logistic regression model showed a small difference between observed and expected probabilities, and there was no significant evidence of poor fit (p-value = 1, greater than 0.05) (Fig. 5J, Table 6 ). Therefore, in logistic regression model with gene-gene interaction, Ho is accepted and the model is accepted as a good fit for the data.

Tumors are commonly treated by surgical removal, chemotherapy, radiotherapy and immunotherapy 57 . Safe surgical removal of glioma is challenging due to its critical location in brain. Additionally, high grade glioma or GBM is highly resistant to chemotherapy and radiotherapy. Immunotherapy treatments are FDA approved for certain blood cancers, but for solid tumors such as GBM immunotherapy ineffective due to incomplete infiltration of immunotherapeutic agent and immune suppression by tumor microenvironment 58, 59 . Diagnosis is the first 1  c1_HuAgeGBsplit_13  1,483  81  c1_HuAgeGBsplit_76  256  161  c1_HuAgeGBsplit_331  143   2  c1_HuAgeGBsplit_26  1,409  82  c1_HuAgeGBsplit_1421  255  162  c1_HuAgeGBsplit_293  141   3  c1_HuAgeGBsplit_10  1,407  83  c1_HuAgeGBsplit_1905  254  163  c1_HuAgeGBsplit_681  141   4  c1_HuAgeGBsplit_6  1,385  84  c1_HuAgeGBsplit_90  242  164  c1_HuAgeGBsplit_1417  140   5  c1_HuAgeGBsplit_180  1,182  85  c1_HuAgeGBsplit_48  241  165  c1_HuAgeGBsplit_212  140   6  c1_HuAgeGBsplit_87  1,144  86  c1_HuAgeGBsplit_1364  234  166  c1_HuAgeGBsplit_84  140   7  c1_HuAgeGBsplit_11  1,070  87  c1_HuAgeGBsplit_294  228  167  c1_HuAgeGBsplit_367  138   8  c1_HuAgeGBsplit_450  969  88  c1_HuAgeGBsplit_58  228  168  c1_HuAgeGBsplit_161  137   9  c1_HuAgeGBsplit_83  930  89  c1_HuAgeGBsplit_179  227  169  c1_HuAgeGBsplit_1195  136   10  c1_HuAgeGBsplit_15  907  90  c1_HuAgeGBsplit_59  224  170  c1_HuAgeGBsplit_1590  136   11  c1_HuAgeGBsplit_22  899  91  c1_HuAgeGBsplit_191  223  171  c1_HuAgeGBsplit_1359  134   12  c1_HuAgeGBsplit_21  875  92  c1_HuAgeGBsplit_27  223  172  c1_HuAgeGBsplit_287 612  104  c1_HuAgeGBsplit_1105  204  184  c1_HuAgeGBsplit_1211  126   25  c1_HuAgeGBsplit_23  602  105  c1_HuAgeGBsplit_446  204  185  c1_HuAgeGBsplit_155  126   26  c1_HuAgeGBsplit_8  582  106  c1_HuAgeGBsplit_49  202  186  c1_HuAgeGBsplit_341  126   27  c1_HuAgeGBsplit_18  571  107  c1_HuAgeGBsplit_292  201  187  c1_HuAgeGBsplit_103  125,   28  c1_HuAgeGBsplit_74  569  108  c1_HuAgeGBsplit_42  201  188  c1_HuAgeGBsplit_159  125   29  c1_HuAgeGBsplit_343  567  109  c1_HuAgeGBsplit_352  198  189  c1_HuAgeGBsplit_368  124   30  c1_HuAgeGBsplit_98  535  110  c1_HuAgeGBsplit_395  198  190  c1_HuAgeGBsplit_4ii  124   31  c1_HuAgeGBsplit_38  528  111  c1_HuAgeGBsplit_156  196  191  c1_HuAgeGBsplit_1614  123   32  c1_HuAgeGBsplit_57  524  112  c1_HuAgeGBsplit_50  194  192  c1_HuAgeGBsplit_459  122   33  c1_HuAgeGBsplit_53  519  113  c1_HuAgeGBsplit_332  192  193  c1_HuAgeGBsplit_776  122   34  c1_HuAgeGBsplit_19  517  114  c1_HuAgeGBsplit_215  190  194  c1_HuAgeGBsplit_1129  121   35  c1_HuAgeGBsplit_188  500  115  c1_HuAgeGBsplit_2103  185  195  c1_HuAgeGBsplit_830  120   36  c1_HuAgeGBsplit_360  487  116  c1_HuAgeGBsplit_2024  184  196  c1_HuAgeGBsplit_107  119   37  c1_HuAgeGBsplit_71  487  117  c1_HuAgeGBsplit_106  183  197  c1_HuAgeGBsplit_164  119   38  c1_HuAgeGBsplit_528  484  118  c1_HuAgeGBsplit_51  182  198  c1_HuAgeGBsplit_1032  117   39  c1_HuAgeGBsplit_29  479  119  c1_HuAgeGBsplit_160  180  199  c1_HuAgeGBsplit_178  117   40  c1_HuAgeGBsplit_103  455  120  c1_HuAgeGBsplit_79  180  200  c1_HuAgeGBsplit_764  116   41  c1_HuAgeGBsplit_14  447  121  c1_HuAgeGBsplit_1106  179  201  c1_HuAgeGBsplit_515  115   42  c1_HuAgeGBsplit_12  443  122  c1_HuAgeGBsplit_291  177  202  c1_HuAgeGBsplit_832  115   43  c1_HuAgeGBsplit_94  423  123  c1_HuAgeGBsplit_75  176  203  c1_HuAgeGBsplit_2184  114   44  c1_HuAgeGBsplit_37  421  124  c1_HuAgeGBsplit_41  175  204  c1_HuAgeGBsplit_456  114   45  c1_HuAgeGBsplit_163  412  125  c1_HuAgeGBsplit_473  170  205  c1_HuAgeGBsplit_97  114   46  c1_HuAgeGBsplit_54  405  126  c1_HuAgeGBsplit_70  170  206  c1_HuAgeGBsplit_356  123   47  c1_HuAgeGBsplit_388  389  127  c1_HuAgeGBsplit_101  167  207  c1_HuAgeGBsplit_295  111   48  c1_HuAgeGBsplit_186  383  128  c1_HuAgeGBsplit_339  167  208  c1_HuAgeGBsplit_470  111   49  c1_HuAgeGBsplit_36  353  129  c1_HuAgeGBsplit_192  166  209  c1_HuAgeGBsplit_757  111   50  c1_HuAgeGBsplit_364  347  130  c1_HuAgeGBsplit_428  166  210  c1_HuAgeGBsplit_153  110   51  c1_HuAgeGBsplit_33  345  131  c1_HuAgeGBsplit_56  164  211  c1_HuAgeGBsplit_549  110   52  c1_HuAgeGBsplit_60  340  132  c1_HuAgeGBsplit_96  163  212  c1_HuAgeGBsplit_44  110   53  c1_HuAgeGBsplit_1177  339  133  c1_HuAgeGBsplit_189  162  213  c1_HuAgeGBsplit_86 www.nature.com/scientificreports/ step in development of an effective treatment plan for any disease, including glioma. GBM is the most aggressive form of glioma that advances quickly giving healthcare providers limited time for diagnosis and treatment 60 . Therefore, to gain deeper understanding of glioma, especially GBM, with hope to develop early accurate diagnosis tools and novel therapies, molecular profiling and network medicine have emerged as research forerunners. Molecular profile or gene based classification and diagnosis help physicians plan treatment and predict clinical outcome. For example, IDH1 gene mutation is highly correlated with glioma survival and is therefore used for glioma classification 61 . IDH1 mutation is lowest in grade 1 glioma that correlates with slow tumor growth and good survival 4 . On the other hand, IDH1 mutation is highest in grade 4 glioma or GBM that correlates with fast tumor growth and poor survival 4 . Improvement in technology, reduced cost of high-throughput sequencing, extensive collaborations and data sharing have made a plethora of glioma molecular profiling datasets such as RNA-seq available to research community 18, 20 . With the explosion of molecular profiling genomics big data, research focus has now shifted from big data mining to big data analysis to prioritize a set of genes with diagnostic value that would eliminate need to profile all 23 K protein coding genes from glioma samples. However, prioritization of a subset of genes for glioma diagnosis and classification has been challenging due to high cellular heterogeneity across and within tumor samples of glioma 62 . Stem cell-like cells in glioma are thought to be responsible for tumor initiation, progression and recurrence 63 . Chemotherapy and radiotherapy kill proliferative stem cells, but are unable to kill quiescent stem cells in the tumor. Quiescent stem cells left in the tumor at end of treatment enter a proliferative state and reconstitute tumor, which leads to tumor recurrence 64 . Therefore, here the goal was to identify distinct sets of proliferative and quiescent stem cell marker genes in GBM that can be used for diagnosis and can serve as potential drug targets.

A meta-analysis was performed on publicly available high-throughput gene expression datasets from human glioma samples, control human brains, normal stem cells and GBM cells in quiescent and proliferative states 18, [26] [27] [28] [29] [30] [31] [32] [33] [34] . DEGs specific to stem cell states were identified from normal stem cell and Glioblastoma cell culture datasets in quiescent and proliferative states. Interestingly, only 45.87% and 25.97% of genes from GBM cell cultures in proliferative and quiescent states were common with normal stem cells in proliferative and quiescent states, respectively (Supplementary Table S1 A,B). This suggests that cancer stem cells, especially those in quiescent state are distinctly different from normal stem cells.

Network analysis facilitates grouping of genes with highly correlated gene expression patterns into modules. It is assumed that modules corelate with distinct biological and cellular states, such as diseases and cell types, respectively. Presently, network analysis identified 9 WGCNA modules and 20 MEGENA modules that were highly correlated with GBM (Fig. 2B,F) . One of these WGCNA modules, HuAgeGBsplit_18 WGCNA module (equivalent c1_HuAgeGBsplit_32/193 MEGENA modules) was also significantly enriched with adult hippocampal rodent quiescent stem cell genes (GSE70696_QNPbyTAP) (Fig. 4B) . Interestingly, though this quiescent stem cell marker enriched HuAgeGBsplit_18 WGCNA module (equivalent c1_HuAgeGBsplit_32/193 MEGENA 56  c1_HuAgeGBsplit_85  334  136  c1_HuAgeGBsplit_162  160  216  c1_HuAgeGBsplit_405  108   57  c1_HuAgeGBsplit_223  323  137  c1_HuAgeGBsplit_217  159  217  c1_HuAgeGBsplit_576  108   58  c1_HuAgeGBsplit_95  321  138  c1_HuAgeGBsplit_2436  159  218  c1_HuAgeGBsplit_1132  106   59  c1_HuAgeGBsplit_501  320  139  c1_HuAgeGBsplit_372  157  219  c1_HuAgeGBsplit_1353  106   60  c1_HuAgeGBsplit_586  308  140  c1_HuAgeGBsplit_477  157  220  c1_HuAgeGBsplit_342  106   61  c1_HuAgeGBsplit_475  307  141  c1_HuAgeGBsplit_68  157  221  c1_HuAgeGBsplit_52  106   62  c1_HuAgeGBsplit_1114  304  142  c1_HuAgeGBsplit_62  156  222  c1_HuAgeGBsplit_1147  105   63  c1_HuAgeGBsplit_1293  303  143  c1_HuAgeGBsplit_151  155  223  c1_HuAgeGBsplit_125  105   64  c1_HuAgeGBsplit_32  302  144  c1_HuAgeGBsplit_622  155  224  c1_HuAgeGBsplit_228  105   65  c1_HuAgeGBsplit_73  298  145  c1_HuAgeGBsplit_1118  152  225  c1_HuAgeGBsplit_361  105   66  c1_HuAgeGBsplit_28  293  146  c1_HuAgeGBsplit_46  152  226  c1_HuAgeGBsplit_553  105   67  c1_HuAgeGBsplit_43  293  147  c1_HuAgeGBsplit_518  151  227  c1_HuAgeGBsplit_2034 www.nature.com/scientificreports/ modules) had a significant correlation with GBM (p-value 0.007) it had a small correlation value of 0.14 with GBM (Fig. 2B) . Possible reasons for this small but significant correlation are discussed here: (A) The result is consistent between WGCNA and MEGENA, two completely different network analysis algorithms. This supports that the results are biologically robust and not a computational artifact that would alter based on alterations in default algorithm settings. (B) SVA + LM normalization was used in this study to retain effects of glioma on gene expression and remove effects of all other covariants. It is possible that effects of covariants such as batch effects are not completely removed by this normalization, which is confounding glioma gene expression effects. (C) It is possible there are other covariants that significantly effect gene expression, such as patients' comorbidities. However, as this information was not available, it could not be included in SVA + LM normalization and therefore glioma effects could not be effectively retained. (D) Controls used in this study comprise of RNA-seq datasets from different parts of the brain derived from humans other than the patients themselves. Ideally control tissue should be derived from the same patient who has the glioma, but presently such patient matched controls were not available for analysis. Lack of patient matched controls is a common challenge in the field of glioma and human disease research. Quiescent stem cells exist in non-proliferative G0 cell cycle phase, but retain ability to reversibly enter cell cycle in response to stimuli. Depletion of surrounding proliferative stem cells and differentiated cells stimulate quiescent stem cells, which are multipotent and have self-renewal potential, to enter cell cycle and replenish the tissue 65 . Quiescent stem cell properties of rodent hippocampal stem cells have been extensively experimentally characterized 33, 66 . Though proliferative stem cell marker genes of GBM origin (PGBCbyQGBC) were significantly enriched in HuAgeGBsplit_18 WGCNA GBM module (equivalent c1_HuAgeGBsplit_32/193 MEGENA GBM modules) with p-value 1.62E-07, it was normal quiescent stem cell marker genes from rodent hippocampal stem cells (GSE70696_QNPbyTAP) that were most significantly enriched with p-value 4.39E-20 (Fig. 4B) . Normal stem cells transition from quiescent state to proliferative state and further to differentiated state. Recently, a set of 336 genes were identified in GBM with SWIM network analysis method, which are potentially involved in transition from stem-like state to differentiated state 52, 53 . Interestingly, no significant enrichment (p-value 0.411, overlap of 9 genes) of SWIM 336 Glioblastoma gene list was found in HuAgeGBsplit_18 GBM WGCNA module (equivalent c1_HuAgeGBsplit_32/193 MEGENA GBM modules). This supports that quiescent stem cell marker genes (GSE70696_QNPbyTAP) enriched in HuAgeGBsplit_18 WGCNA GBM module (equivalent c1_HuAgeGBsplit_32/193 MEGENA GBM modules) represent an undifferentiated quiescent stem cell state distinct from differentiating stem cells.

Gene Ontology (GO) analysis of GSE70696_QNPbyTAP enriched HuAgeGBsplit_18 WGCNA GBM module (equivalent c1_HuAgeGBsplit_32/193 MEGENA GBM modules) revealed enrichment of biomolecule synthesis Table 3 . Common stem cell markers of proliferation.

Common to DEGs listed below from normal stem cells and Glioblastoma in proliferative conditions www.nature.com/scientificreports/ GO terms, such as lipid metabolism, DNA modification, protein post-translational modification, ribosome RNA processing, cell cycle GO terms (G2/M transition, mitotic cell cycle) and signaling pathways such as Wnt signaling ( Fig. 4C-E) . This is consistent with cell cycle and ribosome biogenesis GO terms, previously reported in the rodent hippocampal stem cell dataset from which GSE70696_QNPbyTAP signature genes were identified 33 . As quiescent stem cells can replenish tumor after proliferative stem cells are killed by chemotherapy and radiotherapy, a combinatorial therapy that targets both proliferative and quiescent stem cells could be more effective in GBM treatment. Quiescent stem cell marker genes (GSE70696_QNPbyTAP) enriched in HuAgeGBsplit_18 WGCNA GBM module (equivalent c1_HuAgeGBsplit_32/193 MEGENA GBM modules) reported here, could serve as potential GBM quiescent stem cell drug targets. Small molecule DYRK1B inhibitors were recently shown to target quiescent stem cells and potentiate treatment benefits of chemotherapy 67 . This supports development of treatment strategies to target both proliferative and quiescent stem cells in GBM. Gene expression of eight genes (CD151, CEND1, DCHS1, SMPD1, TPP1, GATD1, RNH1 and SMCR8) from GSE70696_QNPbyTAP enriched in HuAgeGBsplit_18 WGCNA GBM module (equivalent c1_HuAgeG-Bsplit_32/193 MEGENA GBM modules), were sufficient to build a logistic regression diagnostic model that could distinguish between GBM and control samples (Fig. 5J) . Four of the eight genes (CD151, CEND1, SMPD1 and RNH1) used in the model have previously been reported to be important in GBM, other types of cancer or in development [68] [69] [70] [71] . CD151 is a member of tetraspanins scaffolding protein family that is involved in cell-cell adhesion, integrin interaction, cell signaling, cancer progression and metastasis 72 . In GBM, CD151 associates with α3β1 integrin to potentiate EGFR signaling, drive cancer cell motility and tumor aggressiveness 71 . CEND1 or BM88 protein acts as a cell-cycle inhibitor to negatively regulated proliferation and promote differentiation in spinal cord development 70 . In cell membrane lipid bilayer, apoptosis and cellular growth is regulated by balance between sphingosine-1-phosphate and ceramide molecules 73 . SMPD1 gene that regulates 'ceramide sphingosine-1-phosphate rheostat' drives tumor growth and immune escape in non-small cell lung cancer 69 . RNH1 repairs DNA in response to DNA damage and is a known diagnostic and prognostic marker of glioma 68 . Taken together, this provides support for biological relevance of the eight genes (CD151, CEND1, DCHS1, SMPD1, TPP1, GATD1, RNH1 and SMCR8) used here to build a logistic regression diagnostic model for GBM. A molecular Table 5 . Hosmer and Lemeshow goodness of fit (GOF) test for model without gene-gene interaction shows significant difference between observed and expected probabilities (p-value = 0.006132, < 0.05). Model without gene-gene interaction: DiseaseGrade_4 ~ CD151 + CEND1 + DCHS1 + SMPD1 + TPP1 + GATD1 + RNH1 + SMCR8. www.nature.com/scientificreports/ screening kit could be developed with these eight genes for faster and accurate screening of GBM. However, as results of present study were obtained by using computational meta-analysis alone, further experimental validation is required. Overall, this study deconvolutes highly heterogeneous glioma molecular profiles and provides a new perspective to diagnose and develop therapeutic strategies using a small number of quiescent stem cell markers in GBM.

Glioblastoma RNA-seq datasets, SRP027383 and SRP091303, were obtained from NCBI SRA. Other RNAseq datasets for control brains (GSE67333, GSE64810, GSE100297 and GSE53697) and stem cells (GSE68270, GSE70696, GSE99777, GSE93991 and GSE114574), were obtained from NCBI GEO. The author is most grateful to NCBI SRA and NCBI GEO for storing and making these datasets open access. The author also appreciates the organizations and investigators who generously submitted these datasets for sharing to NCBI GEO and NCBI SRA.

Received: 22 December 2019; Accepted: 12 June 2020

Literature review of spinal cord glioblastoma

Survival in glioblastoma: A review on the impact of treatment modalities

Overview of disease and treatment

The 2016 World Health Organization classification of tumors of the central nervous system: A summary

The 2007 WHO classification of tumours of the central nervous system

Hallmarks of glioblastoma: A systematic review. ESMO Open 1, e000144

Current challenges and opportunities in treating glioblastoma

Glioblastoma treatments: An account of recent industrial developments

The cancer stem cell gamble

Secretion-mediated STAT3 activation promotes self-renewal of glioma stem-like cells during hypoxia

Therapeutic strategies targeting cancer stem cells

Emerging targets for glioblastoma stem cell therapy

Cancer stem cell theory: therapeutic implications for nanomedicine

Cancer stem cell quiescence and plasticity as major challenges in cancer therapy

Network medicine: A network-based approach to human disease

A paradigm shift in medicine: A comprehensive review of network-based approaches

Towards precision medicine: The application of omics technologies in asthma management

Comprehensive RNA-seq transcriptomic profiling in the malignant progression of gliomas

Gene expression omnibus: NCBI gene expression and hybridization array data repository

The cancer genome atlas (TCGA): An immeasurable source of knowledge. Contemp

International Cancer Genome Consortium et al. International network of cancer genome projects

Fast R Functions for robust correlations and hierarchical clustering

WGCNA: An R package for weighted correlation network analysis

Multiscale embedded gene co-expression network analysis

International Nucleotide Sequence Database, C. The sequence read archive

RNA sequence analysis of human huntington disease brain reveals an extensive increase in inflammatory and developmental gene expression

Regulatory consequences of neuronal ELAV-like protein binding to coding and non-coding RNAs in human brain

Transcriptomics profiling of Alzheimer's disease reveal neurovascular defects, altered amyloid-beta homeostasis, and deregulated expression of long noncoding RNAs

Cell-specific and region-specific transcriptomics in the multiple sclerosis model: Focus on astrocytes

Quiescence status of glioblastoma stem-like cells involves remodelling of Ca(2+) signalling and mitochondrial shape

Gene signatures of quiescent glioblastoma cells reveal mesenchymal shift and interactions with niche microenvironment

Lysophosphatidic acid receptor is a functional marker of adult hippocampal precursor cells

REST regulation of gene networks in adult neural stem cells

Distinct molecular signatures of quiescent and activated adult neural stem cells reveal specific interactions with their microenvironment

A Microglial signature directing human aging and neurodegeneration-related gene networks

Integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool

Discovering splice junctions with RNA-Seq

HTSeq-A python framework to work with high-throughput sequencing data

A survey of best practices for RNA-seq data analysis

Elegant Graphics for Data Analysis

Visualization of a Correlation Matrix

Capturing heterogeneity in gene expression studies by surrogate variable analysis

The sva package for removing batch effects and other unwanted variation in high-throughput experiments

Removing batch effects and other unwanted noise from sequencing data

limma powers differential expression analyses for RNA-sequencing and microarray studies

Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways

Is my network module preserved and reproducible?

edgeR: A Bioconductor package for differential expression analysis of digital gene expression data

Controlling the false discovery rate in behavior genetics research

Test and visualize gene overlaps

A package for the generation of highly-customizable Venn and Euler diagrams in R

SWIM tool application to expression data of glioblastoma stem-like cell lines, corresponding primary tumors and conventional glioma cell lines

Computational identification of specific genes for glioblastoma stem-like cells identity

Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool

Weighted distributions and estimation of resource selection probability functions

A comparison of goodness-of-fit tests for the logistic regression model

Glioblastoma under Siege: An overview of current therapeutic strategies

Genomics-guided immunotherapy for precision medicine in cancer

Current landscape of immunotherapy in the treatment of solid tumours, with future opportunities and challenges

Glioblastoma multiforme-An overview

Molecular classification of gliomas based on whole genome gene expression: A systematic report of 225 samples from the Chinese Glioma Cooperative Group

Glioblastoma multiforme: A look inside its heterogeneous nature

Therapeutic strategies targeting glioblastoma stem cells. Recent Pat. Anti-Cancer Drug Discov

Mechanisms, hallmarks, and implications of stem cell quiescence

Signaling through BMPR-IA regulates quiescence and long-term activity of neural stem cells in the adult hippocampus

A wake-up call to quiescent cancer cells-Potential use of DYRK1B inhibitors in cancer therapy

Prognostic value of DNA repair genes based on stratification of glioblastomas

Enhanced acid sphingomyelinase activity drives immune evasion and tumor growth in non-small cell lung carcinoma

BM88/CEND1 coordinates cell cycle exit and differentiation of neuronal precursors

CD151-alpha3beta1 integrin complexes are prognostic markers of glioblastoma and cooperate with EGFR to drive tumor cell motility and invasion

CD151 in cancer progression and metastasis: A complex scenario

Ceramide/sphingosine/sphingosine 1-phosphate metabolism on the cell surface and in the extracellular space

This work was conducted at the author's home using personal resources at a) Sarat Ghosh Garden Road, Calcutta, West Bengal, India 700031 and b) Garland Avenue, Downtown Los Angeles, Los Angeles, U.S., 90017 on personal laptops (MacBook Pro and Windows OS) from 2017 to 2019. The author would like to express deepest gratitude to Springer Nature waivers team and Scientific Reports team for the generous discount on publication cost. Other open source and free resources utilized for this work are acknowledged below. 

Author contributions S.M. ideated and managed the work presented in this paper. S.M. also performed experimental design, wrote the manuscript, prepared the figures/tables and submitted the manuscript.

The authors declare no competing interests.

Supplementary information is available for this paper at https ://doi.org/10.1038/s4159 8-020-67753 -5.Correspondence and requests for materials should be addressed to S.M.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creat iveco mmons .org/licen ses/by/4.0/.