Transcriptomic datasets integrated analysis reveals microglia specific biological roles
To explore the specific biological roles of microglia at the glioma condition, we leveraged a series of published RNA sequencing and microarray datasets from brain microglia and macrophage populations isolated from glioma as well as normal samples (see “Materials and methods”). Considering microglia specific biological characterization, we performed two kinds of differential expression analysis, tumor microglia compared to tumor macrophage (MicT/MacT, two datasets) and tumor microglia compared to normal microglia (MicT/MicN, four datasets). We identified the differentially expression genes based on each datasets from MicT/MacT and MicT/MicN (Fig. 1A, see “Materials and methods”). Take the GSE86573 as an example (see Additional file 1: Fig. S1), many known signatures CXCL13, CCL1 and FFAR2 were differentially expressed in MicT/MacT. CXCL13 and CCL1 are all expressed in spinal astrocytes and previous studies have shown that CXCL13 played key roles in some microglia, macrophages, and endothelial cells after CNS infection [14, 15]. For the MicT/MicN results from GSE86573, HAS3, SOX10, AGT and SPAG6 were identified. SOX10, a transcription factor that interacts with Olig2, is important in non-neoplastic oligodendroglial development, and the dys-regulation of mRNA transcripts and protein expression are identified in a wider variety of CNS glial neoplasms [16, 17]. SPAG6 is incorporated into the central apparatus, showing its unambiguously important roles in stabilizing the axoneme [18]. There exists consistent differentially results among multiple datasets both for MicT/MacT and MicT/MicN (Fig. 1A and Additional file 2: Fig. S2). Notably, 9 up-regulated genes were shared by four datasets in MicT/MicN. As shown in Fig. 1A, the functional results in MicT/MacT contained neurotransmitter metabolic process, forebrain neuron differentiation, and positive regulation of odontogenesis (up-regulated) and immune related terms, such as regulation of T cell activation and positive regulation of lymphocyte activation (down-regulated). In MicT/MicN up-regulated results, some developmental and neural terms were involved, including metanephric mesenchyme development (shared by 3 datasets) and regulation of axon guidance (shared by 2 datasets). Neuron differentiation related terms were identified from MicT/MicN down-regulated results. Moreover, some tumor hallmarkers were also related with microglia specific genes, especially the MicT/MicN up-regulation genes (Fig. 1B).
To further explore the relationship among these microglia specific signature genes, we constructed a sub-network based on consistent differential expression genes which were shared by at least two datasets from MicT/MacT and MicT/MicN. Firstly, we obtained a high-quality interaction which were shared by at least two resources based on global protein–protein interaction (PPI) network derived from a previous study [19]. Then, a direct interaction among consistent signature genes based on PPI relation were constructed (Fig. 1C). The interactions with same directions within MicT/MacT or MicT/MicN were defined as “inside the group”, and the other interactions were defined as “outside the group”. As a result, many MicT/MicN signatures were located in the center of this sub-network and closely interacted with other genes, such as JUN, MAP2K4, and LRRK2. And some MicT/MacT signatures displayed “outside the group” interactions, such as PHGDH and MAGI1.
9 MicT/MicN core gene signatures in glioma biology
The 9 common genes shared by four datasets in MicT/MicN were defined as MicT/MicN core genes, and we further dissected the involvement of these microglia specific genes in glioma biological issues, including tumor formation, recurrence and prognosis. A total of 34 public datasets derived from GEO, TCGA, CGGA and PCAWG were obtained (see “Materials and methods”, Additional file 15: Table S1). As shown in Fig. 2A, these 9 MicT/MicN core genes were closely associated with glioma biology. Some genes displayed consistent expression pattern across different glioma events. For example, OAS3 and MMP19 displayed higher expression levels in glioma samples compared to normal samples, recurrence samples compared to primary samples, high-risk group compared to low-risk group (Fig. 2B). MMP19, as a member of the MMPs family, had been reported to promote tumor growth, invasion, metastasis and chemoresistance [20, 21]. And previous study had revealed that glioma patients with higher expression of MMP19 protein had shorter overall survival times [22], which was consistent with our findings. At the aspect of recurrence, the P2RY2 gene displayed higher expression level within three datasets, which was also consistent with the prognostic results.
To test the cell source from immune cells or tumor cells, we further explored the associations between these MicT/MicN core genes and tumor purity, where were estimated by four methods based on TCGA GBM and LGG samples [23]. Among these genes, five genes displayed significant associations with tumor purity (see Fig. 2C and D). Notably, MMP19, LY9 and P2RY2 displayed negative association GBM and LGG purity. Previous researches had found that tumor purity was a prognostic factor for glioma samples [24]. And then, we divided all samples into low-purity and high-purity samples, and the prognostic performance of purity associated genes within these two groups were tested. As shown in Fig. 2D, most genes displayed good predictive performance in LGG than GBM, reflecting the high-grade of GBM. Still, for both GBM and LGG samples, P2RY2 displayed consistent risk predictive performance for both high-purity and low-purity groups, based on four methods (Fig. 2D). Furthermore, we tested the involvement of these P2RY2 signature in clinical response within LGG type with adequate sample number (Fig. 2E). Consistent with the prognostic performance, the risking genes P2RY2 also displayed higher expression level in PD group than SD group (P-value = 0.008), PR group (P-value < 0.001) and CR group (P-value = 0.003). To explore the biological roles, it was observed that the expression level of P2RY2 was associated with wound healing function in GBM type and indel neoantigens and aneuploidy functions in LGG type (Fig. 2F).
Identifying microglia specific subpathways and crosstalk network
Subpathways, regions within the whole pathways, displayed closely association with disease formation and progression, which were revealed by our previous studies [25, 26]. We further developed a novel framework to identify microglia specific subpathways based on both MicT/MacT and MicT/MicN conditions (see Additional file 3: Fig. S3, Materials and Methods). In briefs, we firstly applied network-based random walk algorithm to optimize candidate genes for consistent up-regulated and down-regulated signatures of MicT/MacT and MicT/MicN groups. To validate the robust of random walk results, we obtained another independent data set (GSE29949) to test the results of top 10 genes from MicT/MacT group. As shown in Additional file 4: Fig. S4, nine of ten up-regulated genes displayed significant higher expression level in microglia than other cell types, including monocyte and dendritic cell. And for down-regulated genes ranked by random walk, consistent expression patterns were also observed. Secondly, we obtained subpathway list and calculated the subpathway score by comprehensive considering different dysregulated impacts. And a random strategy (5000 times randomization) was applied to evaluate the significance of all subpathways. As a result, a total of 1/34 up-regulated/down-regulated subpathways from MicT/MacT, and 18/46 up-regulated/down-regulated subpathways from MicT/MicN were identified (see Additional file 5: Fig. S5). Notably, four subpathways were shared by MicT/MacT down-regulated and MicT/MicN down-regulated results, including four subpathways from Regulation of Action Cytoskeleton (Path: 04810) pathway. Finally, based on four types of microglia specific subpathways, we utilized gene overlapping to construct a subpathway crosstalk network. As shown in Additional file 6: Fig. S6, the subpathways with same dysregulated direction displayed closely interactions and most of these subpathways (68/82) were down-regulated. Many subpathways from Axon guidance, Natural killer cell mediated cytotoxicity, and Leukocyte transendothelial migration were the hub subpathways, and four shared subpathways were located in the center of this network.
Microglia subpathways associated with glioma biology
Based on the available glioma data sets used in Fig. 2A, we further explored the associations between microglia specific subpathways and glioma biological events (see “Materials and methods”). As shown in Fig. 3A, all these subpathways displayed two types of expression patterns. Most subpathways from Type I displayed risk expression pattern in glioma formation, recurrence and prognosis, whereas the subpathways from Type II displayed protective pattern. Type I subpathways were majorly derived from Immune system and Cellular community classes, and Type II subpathways were derived from Cell motility, and Development and regeneration classes (see Additional file 7: Fig. S7). At the glioma biology, these subpathways displayed closely associations with tumor formation, especially the prognosis condition. Take the path: 04810_15 (from Type II) as an example, this subpathway displayed higher expression pattern in tumor samples not normal samples, whereas the other subpathways from the same whole pathway (path: 04810) displayed consistent higher expression pattern in normal samples. Next, we calculated the NES score for this subpathway by using the ssGSEA method (see Fig. 3B and C). Similar as the expression pattern from previous result, this subpathway indeed displayed higher functional activity in tumor samples (GSE44971, P-value = 6.3E-05), and higher functional activity in recurrence samples (CGGA693, P-value = 6.8E-08).
A subpathway-based risk signature for glioma prognosis
Our findings revealed that subpathways in the network were closely related with glioma prognosis, especially LGG type (see Fig. 3A). Therefore, we further performed a lasso-based strategy to construct a prognostic model based on 1185 glioma samples of integrated microarray datasets (see “Materials and methods”). As a result, 28 microglia specific subpathways with best predictive performance were identified and defined as SubP28 signature (see “Materials and methods”, Additional file 8: Fig. S8). And the detailed coefficient information of SubP28 signature was provided in Additional file 16: Table S2. For the seven independent testing sets from CGGA, TCGA and PCAWG database, the samples with higher SubP28 score displayed consistent poor survival results than samples with lower score (see Fig. 4A). When considering the aspect of tumor purity, we observed that the SubP28 score of LGG samples was negatively related with tumor purity, however the GBM samples did not (see Additional file 9: Fig. S9A and B). Within the low-purity and high-purity LGG samples, the SubP28 score also predicts the patient survival with significant results (Additional file 9: Fig. S9C and S9D). In the meanwhile, the SubP28 displayed independently predictive performance when considering other clinical factors such as grade, age and IDH1 mutation (see Additional file 9: Fig. S9E).
Furthermore, we quantified our SubP28 signature in data from the Ivy Glioblastoma Atlas Project (IGAP), which performed RNA-seq on microdissections of glioma anatomical structures from hematoxylin and eosin (H&E) staining [Ivy Glioblastoma Atlas Project. http://glioblastoma.alleninstitute.org]. The higher functional activity of SubP28 signatures were enriched in samples from the leading edge of invading glioma and in infiltrating tumor regions. And these subpathway signatures displayed lower activity in cellular tumor and pseudopalisading cells (see Fig. 4B). To test the microglia association of our SubP28 signature, we further obtained several microglia related sets from MsigDB database and calculated the significance of overlapping using the hypergeometric test. As shown in Fig. 4C, the SubP28 signature displayed close association with other microglia (MG) signature.
The glioma samples with different subtypes displayed different prognostic results, and the difference of the SubP28 score between samples with GBM molecular subtypes was also observed (see Fig. 4D). The samples with mesenchymal and neural types exhibited higher SubP28 score than other types. For LGG, we further explored the associations between SubP28 score and IDH1 mutation. The results showed that LGG samples with IDH1 wide type displayed higher SubP28 score. In the meanwhile, the SubP28 displayed significant predictive performance in LGG samples with IDH1 wide type (see Additional file 10: Fig. S10). The correlation analyses revealed that the SubP28 score was positively related with macrophage regulation function, and negatively related with CD4 T cells (see Fig. 4E). Notably, the SubP28 score was also positively related with Macrophage M2 type (a tumor promoting factor), which was consistent with its risk maker in the prognostic results.
SubP28 related with microglia states for both GBM and LGG
To further explore the associations between SubP28 signature and microglia cell state, we calculated the microglia score for TCGA GBM and LGG samples by using the ssGSEA method. And the homeostasis marker (CX3CR1, CSF1R, P2RY12, and TMEM119), M1 marker (IL12B, IL12A, IL23A, TNF, NOS2, and CXCL10) and M2 marker (RETNLB, IL10, ARG1, and MRC1) were utilized for correlation analysis. As shown in Fig. 5A, the SubP28 score of LGG samples were positively related with microglia homeostasis condition, which was not for GBM samples. And for microglia M1 and M2 conditions of GBM samples, the SubP28 score was positively related with microglia activity. It implied the possible mechanism that GBM was majorly consisted of M1/M2 microglia than ones of LGG samples, which was consistent with previous findings. When considering the molecular subtypes and IDH1 mutation, we further observed that M1 and M2 specific associations (see Fig. 5B). It was shown that the SubP28 score displayed positive correlation with M1 markers within LGG samples without IDH1 mutation. For GBM molecular subtypes, the positive correlations with M1 and M2 markers were speciality involved in neural types. Furthermore, we utilized two single-cell RNAseq data sets of glioma to test the SubP28 signature in multiple cells. And it was observed that the SubP28 signature was specific involved in microglia and oligodendrocytes (GSE84465 and GSE89567, see Additional file 11: Fig. S11).
Drug-subpathway network reveals novel treatment strategies
To predict the drug sensitivity in high SubP28 or low SubP28, we systematically evaluated the association between SubP28 score and response sensitivity to anti-neoplastic agents. We obtained the drug response data from three resources, (1) Genomics of Drug Sensitivity in Cancer (GDSC) database for GBM and LGG cell lines [27], (2) Human Glioma Cell Culture (HGCC) cohort, which reported an integrated pharmacogenomic analysis of 100 patient-derived GBM cell cultures treated by 1544 drugs [28], (3) the predicted results from Elastic Net prediction model (LENP), which predict the response sensitivity to antineoplastic compounds for TCGA tumor types, including GBM and LGG [29]. Using the half maximal inhibitory concentration (IC50) value, we calculated the correlation relationship between IC50 of drug/molecules and SubP28 score. Combining the correlation results and drug treatment information, we obtained two candidate drug sets, (1) the drugs displayed higher response sensitivity in cell lines with high SubP28 score, (2) the drugs displayed higher sensitivity with low SubP28 score (see Additional file 12: Fig. S12). Among the drug set i, Lapatinib, a drug for BRCA treatment was identified as sensitive molecule within high SubP28 score group. In the meanwhile, Gemcitabine, Methotrexate, and 5-Fluorouracil were identified to be more sensitive in cases with low SubP28 score.
To explore the detailed associations between antineoplastic compounds and 28 microglia subpathways, we constructed a multi-omic integrated network based on HGCC resource (see “Materials and methods”). As shown in Fig. 6, many microglia specific subpathways were targeted by many candidate drugs from three omics levels, such as path: 03013_19, path: 03008_9, path: 04621_7, and path: 04010_8. And most of these subpathways belonged to the risk subpathway set within the SubP28 model. And the different subpathways (_7 and _8) from the same pathway (path: 04621) displaying opposite prognostic pattern also shared some drugs, such as Clofarabine, Dasatinib, and Thapsigargin. Notably, the path: 04010_8 was enriched by GBM related genes, and GBM mutation genes, including RAC2, NFKB1, RAC1 (disease gene) and MAP4K3, and MAP3K1 (mutation). Also the risk path: 04670_9 was also enriched by GBM related genes, showing its potential roles as the treatment target in the clinical trails.