Quiescent stem cell marker genes in glioma gene networks are sufficient to distinguish between normal and glioblastoma (GBM) samples

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.

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 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

Results
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.

Proliferative and quiescent stem cell marker genes underly GBM modules. Hypergeometric test
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.
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_Symbol
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.

Data availability
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.