- Research article
- Open Access
Subgroup-specific prognostic signaling and metabolic pathways in pediatric medulloblastoma
BMC Cancer volume 19, Article number: 571 (2019)
Using a pathway-focused approach, we aimed to provide a subgroup-specific basis for finding novel therapeutic strategies and further refinement of the risk stratification in pediatric medulloblastoma.
Based on genome-wide Cox regression and Gene Set Enrichment Analysis, we investigated prognosis-related signaling pathways and core genes in pediatric medulloblastoma subgroups using 530 patient data from Medulloblastoma Advanced Genomic International Consortium (MAGIC) project. We further examined the relationship between expression of the prognostic core genes and frequent chromosome aberrations using broad range copy number change data.
In SHH subgroup, relatively high expression of the core genes involved in p53, PLK1, FOXM1, and Aurora B signaling pathways are associated with poor prognosis, and their average expression synergistically increases with co-occurrence of losses of 17p, 14q, or 10q, or gain of 17q. In Group 3, in addition to high MYC expression, relatively elevated expression of PDGFRA, IGF1R, and FGF2 and their downstream genes in PI3K/AKT and MAPK/ERK pathways are related to poor survival outcome, and their average expression is increased with the presence of isochromosome 17q [i(17q)] and synergistically down-regulated with simultaneous losses of 16p, 8q, or 4q. In Group 4, up-regulation of the genes encoding various immune receptors and those involved in NOTCH, NF-κB, PI3K/AKT, or RHOA signaling pathways are associated with worse prognosis. Additionally, the expressions of Notch genes correlate with those of the prognostic immune receptors. Besides the Group 4 patients with previously known prognostic aberration, loss of chromosome 11, those with loss of 8q but without i(17q) show excellent survival outcomes and low average expression of the prognostic core genes whereas those harboring 10q loss, 1q gain, or 12q gain accompanied by i(17q) show bad outcomes. Finally, several metabolic pathways known to be reprogrammed in cancer cells are detected as prognostic pathways including glutamate metabolism in SHH subgroup, pentose phosphate pathway and TCA cycle in Group 3, and folate-mediated one carbon-metabolism in Group 4.
The results underscore several subgroup-specific pathways for potential therapeutic interventions: SHH-GLI-FOXM1 pathway in SHH subgroup, receptor tyrosine kinases and their downstream pathways in Group 3, and immune and inflammatory pathways in Group 4.
With advancements in molecular genomics, many aspects of the genetic basis of tumorigenesis in various cancers have been clarified. Brain tumors have also been a topic of genomic evaluation, and medulloblastoma is one of the most actively studied entities. Consensus has been reached that there are at least 4 subgroups of medulloblastoma: WNT, SHH, Group 3, and Group 4. [1,2,3,4]. In the recent update of the 2016 World Health Organization Classification of Tumors of the Central Nervous System (2016 CNS WHO), the subclassification has been incorporated into the diagnosis of medulloblastoma. Furthermore, preliminary studies have focused on the development of new prognostic biomarkers, novel stratification strategies, or targeted therapies for each subgroup and the elucidation of the subgroup-specific dysregulation of signaling pathways, cytogenetic aberrations, and epigenetic deregulation [2, 5,6,7,8]. However, considerable heterogeneity has remained within the subgroups, among which differences in prognosis are the most clinically important. Combinations of clinical features and several genetic, chromosomal markers have been suggested to classify risk groups within the four subgroups. Recent studies based on integrated multi-omics data including DNA methylation profiles have led to identify refined subtypes within each subgroup of medulloblastoma [9, 10]. The newly defined subtypes were enriched with distinctive clinical and genomic features and several new subtypes were associated with relatively favorable or worse prognosis. Hence, we aimed to discover the biological mechanisms underlying differential clinical outcomes within medulloblastoma subgroups by employing subgroup-specific Gene Set Enrichment Analysis (GSEA) combined with disease progression data analysis . Through the identification of signaling pathways that reflect poor clinical outcome subgroup specifically, we hope to develop new candidates for biomarkers or therapeutic targets for application in precision medicine.
Total 763 expression data of primary medulloblastomas collected via Medulloblastoma Advanced Genomic International Consortium (MAGIC) in a previous study was downloaded from Gene Expression Omnibus (GSE85217) . The raw cel files were preprocessed using the Robust Multi-array Average algorithm  and log2 transformed. To remove multiple probe sets for a given gene in the microarray data, the probe set with the largest inter-quartile range across the samples was selected as a representative one. Sample information including clinical data, subgroup and subtype information, and broad range copy number change data was obtained from the previous study . We further selected only the samples from the patients with age < 18 and with available overall survival information. The final dataset comprised of 530 samples including WNT (n = 49), SHH (n = 121), Group 3 (n = 107), and Group 4 (n = 253). Subsequently, the expression values of each gene were rescaled to relative expression values across the 530 samples ranging from 0 to 1. The 530 samples were divided into exploratory (70% of samples) and validation (30% of samples) datasets, based on randomized sampling stratified by subgroup and 2-year survival status (dead, alive, or censored). Using the exploratory dataset, survival analysis was performed to assess the prognostic impact of individual genes on overall survival in WNT subgroup (n = 34), SHH subgroup (n = 85), Group 3 (n = 74), and Group 4 (n = 177). More specifically, genome-wide univariable Cox regression analysis was employed to calculate proportional hazard ratios (HRs) that indicate the effects of one unit increases in relative expressions of individual genes on risk of death. In this survival analysis step, only the genes with high expression (log2 maximum expression ≥7) and large variation (log2 expression range ≥ 1.5) were analyzed in each subgroup: the final numbers of the genes analyzed were 7832 for WNT, 13,059 for SHH, 12,779 for Group 3, and 13,076 for Group 4. However, in subsequent analyses, WNT subgroup was omitted due to the very low frequency of death event (3 out of 49) which was not suitable for application of survival analysis. To test the fundamental assumption of proportional hazards of the Cox model, we assessed goodness-of-fit for proportional hazards model  and found that the assumption is reasonable for most of the genes investigated (94.9% in SHH subgroup, 96.2% in Group 3, 98.8% in Group 4) (Additional file 1: Table S1). Then, GSEA was performed to detect overrepresented canonical pathways associated with poor prognosis based on curated gene sets of MSigDB (C2) with 1000 non-parametric permutations using pre-ranked genes sorted by the natural log-transformed HRs obtained from the Cox regression analysis . The analysis was confined to total 988 gene sets of size ≤50 to avoid too broad pathways. Initial enriched pathways were detected at a nominal p-value < 0.05 in positive and negative directions respectively, and the genes in the leading-edge subset were considered as core genes in each pathway . Many gene sets were identified redundantly since the gene sets of canonical pathways (C2) in GSEA were collected from multiple databases . Accordingly, we focused on the results from three representative databases, Pathway Interaction Database (PID), Kyoto Encyclopedia of Genes and Genomes (KEGG), and BioCarta (BIOCARTA). Then, collective prognostic power of each pathway was evaluated in both exploratory and validation datasets using Cox regression analysis with average relative expression values of the core genes. To further test the robustness of prognostic significance of each pathway, 1000 bootstrap datasets were generated from entire samples by simple random sampling with replacement. Cox regression model was run on each of the bootstrap samples and bootstrap value for each pathway was calculated by the percentage of bootstrap p-values less than 0.05. We also investigated frequently observed chromosome aberrations that were significantly associated with prognosis in each subgroup. Using broad range copy number change data, we applied univariable Cox regression analysis to detect subgroup-specific prognosis-related gain or loss of each chromosome arm that occurred in each subgroup with a high frequency (≥ 10%). According to the estimated HR, the direction of chromosome aberration is designated “risk” (HR > 1) or “protective” (HR < 1). All data analyses were performed using R statistical software (http://www.R-project.org/).
Ten year Kaplan-Meier plots based on 530 overall survival data of pediatric medulloblastoma patients revealed that survival prognoses of the four subgroups were largely congruent with previous studies in both exploratory and validation datasets [2, 5, 8]. The WNT subgroup was associated with the best survival, SHH subgroup and Group 4 with intermediate survival, and Group 3 with the poorest survival (Additional file 2: Figure S1).
In GSEA analysis, nine canonical pathways in SHH subgroup are detected to be prognostic toward positive direction including four PID canonical pathways, “p53 pathway”, “PLK1 signaling events”, “Aurora B signaling”, and “ FOXM1 transcription factor network”, three KEGG pathways, “Alanine, aspartate, and glutamate metabolism”, “RNA polymerase”, and “Basal transcription factors”, and two BIOCARTA pathways, “Cyclins and Cell Cycle Regulation” and “Proteasome complex” (Table 1, Additional file 1: Table S2). In bootstrap analysis with 1000 replications, eight out of nine canonical pathways show statistically significant results in more than 90% of resampled datasets. According to the average expression values of 31 core genes from the top two PID canonical pathways (“p53 pathway” and “PLK1 signaling events”), the patients are separated into three different survival groups similarly in both exploratory and validation datasets (Fig. 1a). Furthermore, we also obtain statically significant results with a whole dataset (n = 121) in log-rank tests with average expression of core genes in nine canonical pathways respectively (Additional file 1: Table S2).
In association analyses between survival outcomes and copy number changes, gain of 9p is detected as a significant protective aberration (Additional file 1: Table S3). The patients with 9p gain show excellent prognosis (Fig. 1b), however, the expression of 31 prognostic core genes is not significantly low in the tumors with 9p gain (Fig. 1b boxplot). On the other hand, losses of 17p, 14q, 10q, and isochromosome 17q [i(17q)] are observed to be risk aberrations (Additional file 1: Table S3). After excluding the patients with 9p gain, we further find that 16 patients harboring i(17q) or more than two losses out of 10q, 14q, and 17p show the worst survival outcome (Fig. 1c). Seven patients with either one of 10q loss or 14q loss show intermediate survival rates while three patients with 17p loss alone show the best prognosis. Furthermore, the highest expression level of 31 prognostic core genes is observed in the 16 patients with i(17q) or multiples losses of 10q, 14q, or 17p whereas the lowest is detected in the three patients with 17p loss alone (Fig. 1c boxplot). Finally, 69 patients who do not possess any risk or protective aberrations can be divided into two different survival groups according to the average expression of the 31 core genes (Fig. 1d). The average expressions of the core genes from the nine prognostic canonical pathways is presented in a heatmap including clinical information, subtype information recently reported by Cavalli et al. , and the stratification of patients according to the status of prognostic chromosome aberrations and gene expressions shown in Fig. 1b-d (Fig. 1e). The heatmap shows that the average expressions of the core genes of the nine pathways considerably correlate with one another. On the other hand, the average expression of the 31 core genes from the top two PID canonical pathways is also significantly associated with subtype, age group, and histology; the highest expression is observed in subtype of SHH α and histological group of large cell/anaplastic (LCA) (Additional file 2: Figure S2A).
Total 35 canonical pathways in PID, KEGG, and BIOCARTA are associated with poor prognosis in Group 3 (Additional file 1: Table S2). Of the 35 canonical pathways, 16 pathways show > 90% of statistical significance in 1000 bootstrap samples, including many cancer-related signaling pathways such as “Notch-mediated HES/HEY network” (Notch network), “Class I PI3K signaling events mediated by Akt” (PI3K/Akt), “p53 pathway”, “Melanoma”, “Bladder cancer”, “Proteasome”, and “Influence of Ras and Rho proteins on G1 to S Transition” (Ras/Rho on G1 to S). In addition, the prognostic canonical pathways include several metabolic pathways such as “Tyrosine and galactose metabolism”, “Pentose phosphate pathway”, and “TCA cycle”. Top 5 prognostic canonical pathways in each database ranked by bootstrap percentage are presented with representative prognostic core genes in Table 2. According to the average expression of 20 core genes from the top two PID canonical pathways (“Notch network” and “PI3K/Akt”), low and high expression groups of the Group 3 patients are clearly separated into different survival groups even though medium groups are inconsistently segregated into the low or high group in exploratory and validation datasets respectively (Fig. 2a).
Several chromosome aberrations including losses of 16p, 8q, 16q, 4q, 22q, and 11q are detected as protective aberrations (Additional file 1: Table S4). However, the protective aberrations are significantly associated with one another and many patients possess multiple protective aberrations simultaneously. Thus, we applied multivariable Cox regression analysis and selected the losses of 16p, 8q and 4q as the final protective aberrations because the losses of 16q, 22q, and 11q contributed little to increase in R-squared value. When 51 patients possessing one of the protective aberrations, losses of 16p, 8q and 4q, are divided into three groups according to the number of the losses, seven patients harboring all the losses show the most favorable prognosis and the survival rates are reduced as the number of losses decrease (Fig. 2b). Besides, the number of losses inversely correlates with the average expression levels of the 20 prognostic core genes from the top two PID pathways (Fig. 2b, boxplot). On the other hand, 17p loss and i(17q) are detected as risk aberrations (Additional file 1: Table S4). In multivariable Cox regression analysis, 17p loss is no longer significant when the status of i(17q) is added to the model. Moreover, the patients with 17p loss without 17q gain do not show difference in prognosis compared to those without 17p loss (p-value = 0.49). In contrast, loss of 17p accompanied by loss of 10q reduces survival rates considerably (Fig. 2c). Hence, we conclude that i(17q) and simultaneous loss of 17p and 10q are the risk aberrations in Group 3. After excluding 77 patients with the protective or risk aberrations, the prognosis of the remained 30 patients is separated according to the expression level of 20 core genes (Fig. 2d). The heatmap of core gene expressions of 15 prognostic gene sets shows that protective and risk chromosome aberrations are associated with low and high average expressions of the core genes respectively (Fig. 2e). The expression of the core genes is also significantly associated with subtype in Group 3; the highest expression is observed in subtype of Group 3 γ and histological group of LCA (Additional file 2: Figure S2B).
In Group 4, total 41 canonical pathways out of 58 canonical pathways significantly detected in GSEA show statistical significance in more than 90% of 1000 bootstrap samples, including “Notch network”, “p53 pathway”, “RhoA signaling pathway”, “NF-κB Signaling Pathway”, “One carbon pool by folate”, “Cell Cycle: G1/S Check Point”, “Integrin-linked kinase signaling”, and several immune response-related canonical pathways such as “TNF receptor signaling pathway”, “Fc-epsilon receptor I signaling in mast cells”, and “Antigen processing and presentation” (Additional file 1: Table S2, Table 3). Stratification of the Group 4 patients using average expression of 31 core genes from the top two PID canonical pathways (“Notch network” and “TNF receptor signaling pathway”) produces clear separation of survival curves in both exploratory and validation datasets (Fig. 3a).
Investigation of prognostic chromosome aberration reveals that losses of 11p and 8q are protective aberrations (Additional file 1: Table S5). In previous studies, loss of chr11 has been recognized as an important cytogenetic marker of favorable prognosis in Group 4 [2, 5, 15]. Congruently, we find that 11p loss not accompanied by 11q loss is not a significant prognostic factor (data not shown), and simultaneous loss of 11p and 11q is the most significant prognostic chromosome aberration in Group 4 (Additional file 1: Table S5). The patients with loss of chr11 show excellent survival outcomes and low average expression of 31 prognostic core genes from the top two PID pathways (Fig. 3b). In multivariable Cox regression analysis, the other protective aberration, 8q loss, is no longer significant (p-value = 0.23) when the status of chr11 is added to the model. Therefore, we conclude that chr11 loss is the strongest protective chromosome aberration in Group 4. We further investigate prognostic chromosome aberrations after excluding the patients with loss of chr11 in Group 4, and find that i(17q) is a risk chromosome aberration (Additional file 1: Table S6, Fig. 3c). In a similar successive manner, we find that 8q loss is a significant protective aberration in 90 patients after excluding 163 Group 4 patients with chr11 loss or i(17q) (Additional file 1: Table S7, Fig. 3d). Likewise, 10q loss, 1q gain, and 12q gain, are identified as risk aberrations within 120 patients with i(17q) but without loss of chr11 (Additional file 1: Table S8, Fig. 3e). Finally, two groups of remained patients, 91 patients with i(17q) and 58 patients without i(17q), can be further separated into different survival groups respectively according to the average expression of 31 prognostic core genes (Fig. 3f and g), which can be also shown by the heatmap with detailed information (Fig. 3h). On the other hand, high expression level of the 31 core genes is observed in the patients with leptomeningeal metastases (M+) (Additional file 2: Figure S2C).
Current consensus of risk stratification of childhood medulloblastoma has refined the risk stratification according to subgroup and an integration of the clinical markers as well as the molecular, genomic profiles of the tumor . In the SHH subgroup, the risk group was divided by the combination of TP53 mutation, MYCN amplification, and metastasis. For Group 3, the major factors deciding risk were metastasis and MYC amplification. In Group 4, the important segregation points were metastasis and chr11 loss. Despite this huge improvement and modification, ‘unknown’ types of patients have remained, such as Group 3 patients who are non-metastatic but have MYC amplification. Furthermore, metastasis is even more prominent as the major factor in all subgroups, and molecular and genomic mechanisms that are masked by the clinical phenotype of metastasis have remained to be elucidated. More recent studies have examined the presence of subtypes in each medulloblastoma subgroup, characterized by different clinical, genomic, and epigenomic features [9, 10, 16], indicating considerable within-subgroup heterogeneity. Thus, we investigate prognostic signaling pathways in each subgroup of medulloblastoma, to provide molecular basis for characterization of within-subgroup heterogeneity and subgroup-specific therapeutic strategies. Finally, the results are comprehensively summarized in schematic illustrations of prognosis-related signaling pathways (Fig. 4) and metabolic pathways (Fig. 5).
In SHH subgroup, majority of the prognostic gene sets are related to cell cycle and mitosis such as the canonical pathways of “p53 pathway”, “PLK1 signaling events”, “FOXM1 transcription factor network”, “Aurora B signaling”, and “Cyclins/CDKs”. In SHH signaling pathway, FOXM1 has been known to be a downstream target gene of GLI transcription factor that is activated by SMO relieved by hedgehog ligands (Fig. 4a) [17,18,19]. Consistently, we observe a considerable correlation between the expressions of GLI2 and FOXM1 (Pearson’s correlation coefficient = 0.48, p-value = 2.9 × 10− 8). Additionally, as expected, the expressions of GLI1 and GLI2 are profoundly high in SHH subgroup (Additional file 2: Figure S3A) and high expression level of GLI2 is significantly associated with prognosis (HR = 1.45 per 0.1 increment of relative expression, p-value = 0.0074). As a transcription factor, FOXM1 subsequently activates essential genes for mitotic progression including PLK1 [20,21,22]. FOXM1 and PLK1 are cooperatively overexpressed in various cancers [23,24,25] and identified as potential therapeutic targets [22, 26,27,28]. They play crucial roles in mitosis entry by complicated inter-activating relationship . Such a reciprocal activation is also reported between PLK1 and MYCN , one of the well-known genes highly expressed in SHH medulloblastomas. Besides the core genes detected in GSEA, two genes encoding ErbB family proteins, EGFR and ERBB3, are observed to be prognostic in gene-wise univariable Cox regressions (HR = 1.35 per 0.1 increment of relative expression and p-value = 0.045 for EGFR; HR = 1.55 per 0.1 increment of relative expression and p-value < 0.00001 for ERBB3). The expression of EGFR is relatively high in SHH subgroup. In contrast, only a small subset of tumors shows a very high level of ERBB3 expression (Additional file 2: Figure S3A).
Gain of 9p, frequently co-occurred with loss of 9q, was recognized in SHH tumors , however, its association with prognosis has not been reported. In our result, favorable prognosis is observed in the SHH patients with gain of 9p. It has been well known that 9p contains two important tumor suppressor genes, CDKN2A and CDKN2B. Indeed, median expression of CDKN2B was slightly higher in SHH patients with 9p gain with statistical significance (Additional file 2: Figure S3B). We also detect that the SHH patients with i(17q) or simultaneous loss of 10q, 17p, or 14q have the worst survival outcome (median survival time = 1.59 years) (Fig. 1c). As expected, the expressions of important tumor suppressors, TP53 and PTEN/SUFU, are significantly or marginally significantly reduced with the corresponding chromosome aberrations, 17p loss and 10q loss in SHH medulloblastomas (Additional file 2: Figure S3B). In Fig. 4a, several anti-cancer agents targeting deregulated signaling pathways in SHH subgroup are presented. Proteasome inhibitors (bortezomib) and thiazole antibiotics (sisomicin A and thiostrepton) target FOXM1 and induce apoptosis [32, 33]. The mitotic kinases PLK1 and AURKA/B are blocked by the small-molecule selective inhibitors BI 2536/volasertib and alisertib, respectively [34, 35]. Palbociclib and abemaciclib are selective CDK4/6 inhibitors . SHH signaling is also blocked by the SMO inhibitors vismodegib and saridegib and the GLI inhibitors GANT58, GANT61, and arsenic trioxide. Up-regulated EGFR is blocked by various inhibitors including gefitinib, erlotinib, and cetuximab.
In Group 3, two receptor tyrosine kinases (PDGFRA and IGF1R), FGF2 also known as basic fibroblast growth factor (bFGF), and the genes in their downstream signaling pathways such as PI3K/AKT and MAPK/ERK are detected as prognostic core genes (Fig. 5a). Furthermore, IGF1 is observed to be most highly expressed in Group 3 (Additional file 2: Figure S4A). High expression of MYC is confirmed in Group 3 tumors, and its target genes involved in glycolysis including SLC2A1 (GLUT1), LDHA, and HK2 are noticed to be highly expressed as well (Additional file 2: Figure S4A). Additionally, many genes encoding proteasomes, cyclins, and CDKs are also detected as prognostic genes.
The most favorable survival rate is observed in the Group 3 patients who have simultaneous loss of 4q, 8q, or 16p, which can be possibly explained by the loss of the oncogenic core genes such as MYC on 8q, FGF2 and PDGFRA on 4q, and PLK1 on 16p. Indeed, the expressions of the four genes are significantly or marginally significantly low in accordance with the corresponding chromosomal deletions (Additional file 2: Figure S4B). On the other hand, i(17q) and simultaneous loss of 17p and 10q are associated with poor survival outcome in Group 3, which is also observed in SHH subgroup.
Anti-cancer agents targeting PDGFRA include receptor tyrosine kinase (RTK) inhibitors such as imatinib, sunitinib, and dasatinib. Downstream of RTKs is targeted by PI3K inhibitors (buparlisib and BYL719) [37, 38] and Ras-Raf-MEK-ERK pathway inhibitors (diazepinomicin, sorafenib and regorafenib) [39,40,41]. Up-regulated proteasome complex is blocked by bortezomib.
In Group 4, “Notch network” genes such as NOTCH1, NOTCH2, and HEY1 are included in the prognostic core genes (Additional file 1: Table S2, Table 3). In addition, high expressions of many genes encoding immune receptors are also associated with poor survival outcome, including cytokine receptors (TNFRSF1A/B, ITGB1/2, IFNAR2, IL1R1, IL4R, IL7R, IL10RA, CXCR4, and CD74), Toll-like receptors (TLR3/4/5), and Fc receptors (FCER1G and FCGR2A/B). Interestingly, the expressions of NOTCH1, NOTCH2 or NOTCH3 are highly correlated with those of many prognostic immune receptors (Additional file 2: Figure S5), implying that Notch signaling might be important in the regulation of tumor immune response or tumor microenvironment, which, in fact, has been reported in many studies [42, 43]. For example, activation of Notch signaling induced recruitment of tumor-associated macrophage (TAM)  and anti-Notch treatment reduced immune-suppressive M2 TAM infiltration  or inhibited tumor-induced T-cell tolerance mediated by Myeloid-derived suppressor cells , thus underscoring possible immunotherapeutic opportunities in Group 4 medulloblastoma patients with relative activation of Notch signaling. Notch signaling is also known to be closely related to stemness of medulloblastoma cell and essential for maintenance of stemness in brain tumor initiating cell [46, 47]. Enhanced expressions of the genes involved in nucleotide synthesis (DHFR and TYMS), cell cycle, and NF-κB, PI3K/AKT, or RHOA pathways are also associated with poor prognosis. Several genes in JAK-STAT pathway, the major downstream signaling of cytokine receptors, are most highly expressed in Group 4, including JAK2, STAT2, STAT3, STAT5A, and STAT5B (Additional file 2: Figure S6A). Besides, “Validated nuclear estrogen receptor α network” is recognized as the deregulated signaling pathway in Group 4 (Table 3). Congruently, expression of estrogen-related receptor γ (ESRRG) is found to be considerably high in Group 4 tumors (Additional file 2: Figure S6A).
The Group 4 patients with loss of chr11 or loss of 8q without i(17q) show excellent survival rates (Fig. 3b and d), which could possibly be explained by reduced expressions of prognostic core genes such as HRAS on 11p, GAB2 and RELA on 11q, and MYC and HEY1 on 8q (Additional file 2: Figure S6B). Likewise, the poor survival outcome of the Group 4 patients with loss of 10q, or gains of 1q or 12q accompanied by i(17q) might be related to reduced expression of PTEN on 10q, or elevated expressions of AKT3 and PSEN2 on 1q and CDK2 on 12q (Additional file 2: Figure S6B, Fig. 4c).
The Notch signaling pathway is blocked by MK-0752 and RO4929097, inhibitors of γ-secretase that mediates the cleavage of the Notch intracellular domain (NICD), which is translocated to the nucleus and activates the transcription of NOTCH target genes [48, 49]. Considering the prognostic significance of immune receptors, NF-κB signaling, and JAK-STAT pathway, anti-inflammatory agents and monoclonal antibodies against immune receptors (milatuzumab and infliximab) may provide therapeutic interventions in Group 4 [50,51,52]. Anti-cancer agents targeting nucleotide biosynthesis pathways include a DHFR inhibitor (methotrexate) and a TYMS inhibitor (5-FU).
Prognostic metabolic pathways in medulloblastoma subgroups
Finally, we also detect several metabolic pathways that are known to be reprogrammed in cancer as prognostic canonical pathways such as “Alanine, aspartate, and glutamate metabolism” in SHH subgroup, “Pentose phosphate pathway” and “TCA cycle” in Group 3, and “One carbon pool by folate” in Group 4. Accordingly, focusing on the widely known metabolic pathways modified in cancer cells [53,54,55,56,57], a comprehensive schematic diagram is presented in Fig. 5 based on not only the prognostic core genes obtained in GSEA but also the genes significantly detected in gene-wise univariable Cox regression analysis (p-value < 0.05) or highly expressed in a particular medulloblastoma subgroup. Several genes are commonly detected as prognostic genes in all three subgroups, SHH subgroup, Group 3, and Group 4, including ENO1 encoding enolase 1 in glycolytic pathway and FASN and SCD encoding fatty acid synthase and stearoyl-CoA desaturase respectively in fatty acid synthesis pathway. Many genes involved in glycolysis or glutaminolysis are prognostic or highly expressed in SHH subgroup or Group 3. The genes involved in one carbon (folate) cycle or serine synthesis are prognostic mainly in SHH subgroup or Group 4. Two important genes in pentose phosphate pathway, TKT and TALDO1, are associated with poor prognosis in Group 3 and Group 4. It has long been recognized that activation of those metabolic pathways supports tumor cell anabolism for rapid proliferation and plays a crucial role in maintenance of energy and redox homeostasis in cancer cells [53,54,55,56,57,58]. Thus, targeting the key enzymes in the altered tumor metabolic pathways is now being actively investigated as a new potential chemotherapeutic strategy [59, 60], which also shed light on treatment for medulloblastoma.
Data from a decade of genomic analysis have completely changed the definition and diagnosis of medulloblastoma. Research is now progressing toward the identification of prognostic and therapeutic applications of this accumulated genomic data. Our primary and ultimate goal is to identify the genes or signaling pathways that are essentially responsible for the prognosis of each established medulloblastoma subgroup. The prognostic genes detected in a particular subgroup are not necessarily highly expressed in the given subgroup compared to other subgroups. Rather, we find that most of the subgroup-specific prognosis-related genes show almost the same levels or even lower levels of expression, indicating complicated prognostic molecular determinants. Thus, our results suggest that we need to consider not only the genes that are highly expressed or amplified but also those with average or relatively low level of expression to expand the therapeutic targets and strategies.
Availability of data and materials
The datasets analyzed during the current study are available in the Gene Expression Omnibus (GSE85217).
Acetyl-CoA Carboxylase Alpha
ATP Citrate Lyase
AKT Serine/Threonine Kinase
Aldolase Fructose-Bisphosphate A
Aurora kinase A
basic fibroblast growth factor
Baculoviral IAP Repeat Containing 5 (Survivin)
Centromere protein F
C-X-C Motif Chemokine Receptor 4
E2F Transcription Factor
Epidermal Growth Factor Receptor
Erb-B2 Receptor Tyrosine Kinase 3
Extracellular signal–regulated kinase
Fatty Acid Synthase
Fibroblast growth factor 2
Fibroblast growth factor receptor
Fos proto-oncogene AP-1 transcription factor subunit
Forkhead box protein M1
GRB2 Associated Binding Protein 2
phosphoribosylglycinamide formyltransferase phosphoribosylglycinamide synthetase phosphoribosylaminoimidazole synthetase
GLI Family Zinc Finger
Glutamate Dehydrogenase 1
Glutamic-Oxaloacetic Transaminase 2
Glutamic--Pyruvic Transaminase 2
Gene Set Enrichment Analysis
Hes Related Family BHLH Transcription Factor With YRPW Motif 1
3-Hydroxy-3-Methylglutaryl-CoA Synthase 1
HRas Proto-Oncogene GTPase
Isocitrate Dehydrogenase (NADP(+)) 1 Cytosolic
Isocitrate Dehydrogenase (NADP(+)) 2 Mitochondrial
Interferon Alpha And Beta Receptor Subunit 2
Insulin Like Growth Factor 1
Insulin Like Growth Factor 1 Receptor
NFKB Inhibitor Alpha (NFKBIA)
Interleukin 1 Receptor Type 1
Interleukin 4 Receptor
Interleukin 7 Receptor
Kyoto Encyclopedia of Genes and Genomes
KRAS Proto-Oncogene GTPase
Aldehyde Dehydrogenase 1 Family Member L1
Lactate Dehydrogenase A
Medulloblastoma Advanced Genomic International Consortium
Mitogen-Activated Protein Kinase
Malate Dehydrogenase 1
Malate Dehydrogenase 2
Methylenetetrahydrofolate Dehydrogenase (NADP+ Dependent) 1 Like
Methylenetetrahydrofolate Dehydrogenase (NADP+ Dependent) 2 Methenyltetrahydrofolate Cyclohydrolase
Methylenetetrahydrofolate Dehydrogenase (NADP+ Dependent) 2 Like
MYC Proto-Oncogene BHLH Transcription Factor
MYCN Proto-Oncogene BHLH Transcription Factor
NIMA Related Kinase 2
Notch intracellular domain
Nuclear Factor Kappa B Subunit 1 (NFKB1)
Platelet Derived Growth Factor Receptor Alpha
3-Phosphoinositide Dependent Protein Kinase 1
Phosphoglycerate Mutase 2
Class I PI3K signaling events mediated by Akt
Pathway Interaction Database
Polo-like kinase 1
Phosphoribosyl Pyrophosphate Amidotransferase
Phosphoribosyl Pyrophosphate Synthetase 2
Phosphoserine Aminotransferase 1
Raf-1 Proto-Oncogene Serine/Threonine Kinase
- Ras/Rho on G1 to S:
Influence of Ras and Rho proteins on G1 to S Transition
RELA Proto-Oncogene NF-KB Subunit
Ras Homolog Family Member A
Ribosomal Protein S6 Kinase B1 (p70S6K)
Succinate Dehydrogenase Complex Iron Sulfur Subunit B
Serine Hydroxymethyltransferase 2
Solute Carrier Family 16 Member 1
Solute Carrier Family 1 Member 5
Solute Carrier Family 2 Member 1 (GLUT1)
Solute Carrier Family 2 Member 3 (GLUT3)
Solute Carrier Family 38 Member 5
Smoothened Frizzled Class Receptor
signal transducer and activator of transcription
Succinate-CoA Ligase Alpha Subunit
SUFU Negative Regulator of Hedgehog Signaling
Toll-like receptor 4
tumor necrosis factor receptor
Northcott PA, Dubuc AM, Pfister S, Taylor MD. Molecular subgroups of medulloblastoma. Expert Rev Neurother. 2012;12(7):871–84.
Kool M, Korshunov A, Remke M, Jones DT, Schlanstein M, Northcott PA, Cho YJ, Koster J, Schouten-van Meeteren A, van Vuurden D, et al. Molecular subgroups of medulloblastoma: an international meta-analysis of transcriptome, genetic aberrations, and clinical data of WNT, SHH, group 3, and group 4 medulloblastomas. Acta Neuropathol. 2012;123(4):473–84.
Taylor MD, Northcott PA, Korshunov A, Remke M, Cho YJ, Clifford SC, Eberhart CG, Parsons DW, Rutkowski S, Gajjar A, et al. Molecular subgroups of medulloblastoma: the current consensus. Acta Neuropathol. 2012;123(4):465–72.
Li KK, Lau KM, Ng HK. Signaling pathway and molecular subgroups of medulloblastoma. Int J Clin Exp Pathol. 2013;6(7):1211–22.
Shih DJ, Northcott PA, Remke M, Korshunov A, Ramaswamy V, Kool M, Luu B, Yao Y, Wang X, Dubuc AM, et al. Cytogenetic prognostication within medulloblastoma subgroups. J Clin Oncol. 2014;32(9):886–96.
MacDonald TJ, Aguilera D, Castellino RC. The rationale for targeted therapies in medulloblastoma. Neuro-Oncology. 2014;16(1):9–20.
Remke M, Ramaswamy V, Taylor MD. Medulloblastoma molecular dissection: the way toward targeted therapy. Curr Opin Oncol. 2013;25(6):674–81.
Northcott PA, Jones DT, Kool M, Robinson GW, Gilbertson RJ, Cho YJ, Pomeroy SL, Korshunov A, Lichter P, Taylor MD, et al. Medulloblastomics: the end of the beginning. Nat Rev Cancer. 2012;12(12):818–34.
Cavalli FMG, Remke M, Rampasek L, Peacock J, Shih DJH, Luu B, Garzia L, Torchia J, Nor C, Morrissy AS, et al. Intertumoral heterogeneity within Medulloblastoma subgroups. Cancer Cell. 2017;31(6):737–754 e736.
Schwalbe EC, Lindsey JC, Nakjang S, Crosier S, Smith AJ, Hicks D, Rafiee G, Hill RM, Iliasova A, Stone T, et al. Novel molecular subgroups for clinical classification and outcome prediction in childhood medulloblastoma: a cohort study. Lancet Oncol. 2017;18(7):958–71.
Park AK, Kim P, Ballester LY, Esquenazi Y, Zhao Z. Subtype-specific signaling pathways and genomic aberrations associated with prognosis of glioblastoma. Neuro Oncol. 2019;21(1):59-70.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4(2):249–64.
Chang CC, Weissfeld LA. Normal approximation diagnostics for the Cox model. Biometrics. 1999;55(4):1114–9.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Ramaswamy V, Remke M, Bouffet E, Bailey S, Clifford SC, Doz F, Kool M, Dufour C, Vassal G, Milde T, et al. Risk stratification of childhood medulloblastoma in the molecular era: the current consensus. Acta Neuropathol. 2016;131(6):821–31.
Northcott PA, Buchhalter I, Morrissy AS, Hovestadt V, Weischenfeldt J, Ehrenberger T, Grobner S, Segura-Wang M, Zichner T, Rudneva VA, et al. The whole-genome landscape of medulloblastoma subtypes. Nature. 2017;547(7663):311–7.
Quan M, Wang P, Cui J, Gao Y, Xie K. The roles of FOXM1 in pancreatic stem cells and carcinogenesis. Mol Cancer. 2013;12:159.
Teh MT, Wong ST, Neill GW, Ghali LR, Philpott MP, Quinn AG. FOXM1 is a downstream target of Gli1 in basal cell carcinomas. Cancer Res. 2002;62(16):4773–80.
Douard R, Moutereau S, Pernet P, Chimingqi M, Allory Y, Manivet P, Conti M, Vaubourdolle M, Cugnenc PH, Loric S. Sonic hedgehog-dependent proliferation in a series of patients with colorectal cancer. Surgery. 2006;139(5):665–70.
Zhu H. Targeting forkhead box transcription factors FOXM1 and FOXO in leukemia (review). Oncol Rep. 2014;32(4):1327–34.
Costa RH. FoxM1 dances with mitosis. Nat Cell Biol. 2005;7(2):108–10.
Koo CY, Muir KW, Lam EW. FOXM1: from cancer initiation to progression and treatment. Biochim Biophys Acta. 2012;1819(1):28–37.
Dibb M, Han N, Choudhury J, Hayes S, Valentine H, West C, Ang YS, Sharrocks AD. The FOXM1-PLK1 axis is commonly upregulated in oesophageal adenocarcinoma. Br J Cancer. 2012;107(10):1766–75.
Dibb M, Han N, Choudhury J, Hayes S, Valentine H, West C, Sharrocks AD, Ang YS. FOXM1 and polo-like kinase 1 are co-ordinately overexpressed in patients with gastric adenocarcinomas. BMC Res Notes. 2015;8:676.
Zhang Z, Zhang G, Kong C. FOXM1 participates in PLK1-regulated cell cycle progression in renal cell cancer cells. Oncol Lett. 2016;11(4):2685–91.
Wierstra I. FOXM1 (Forkhead box M1) in tumorigenesis: overexpression in human cancer, implication in tumorigenesis, oncogenic functions, tumor-suppressive properties, and target of anticancer therapy. Adv Cancer Res. 2013;119:191–419.
Degenhardt Y, Lampkin T. Targeting polo-like kinase in cancer therapy. Clin Cancer Res. 2010;16(2):384–9.
Strebhardt K, Ullrich A. Targeting polo-like kinase 1 for cancer therapy. Nat Rev Cancer. 2006;6(4):321–30.
Fu Z, Malureanu L, Huang J, Wang W, Li H, van Deursen JM, Tindall DJ, Chen J. Plk1-dependent phosphorylation of FoxM1 regulates a transcriptional programme required for mitotic progression. Nat Cell Biol. 2008;10(9):1076–82.
Xiao D, Yue M, Su H, Ren P, Jiang J, Li F, Hu Y, Du H, Liu H, Qing G. Polo-like Kinase-1 regulates Myc stabilization and activates a feedforward circuit promoting tumor cell survival. Mol Cell. 2016;64(3):493–506.
Northcott PA, Korshunov A, Witt H, Hielscher T, Eberhart CG, Mack S, Bouffet E, Clifford SC, Hawkins CE, French P, et al. Medulloblastoma comprises four distinct molecular variants. J Clin Oncol. 2011;29(11):1408–14.
Bhat UG, Halasi M, Gartel AL. FoxM1 is a general target for proteasome inhibitors. PLoS One. 2009;4(8):e6593.
Liu LL, Zhang DH, Mao X, Zhang XH, Zhang B. Over-expression of FoxM1 is associated with adverse prognosis and FLT3-ITD in acute myeloid leukemia. Biochem Biophys Res Commun. 2014;446(1):280–5.
Gjertsen BT, Schoffski P. Discovery and development of the polo-like kinase inhibitor volasertib in cancer therapy. Leukemia. 2015;29(1):11–9.
Friedberg JW, Mahadevan D, Cebula E, Persky D, Lossos I, Agarwal AB, Jung J, Burack R, Zhou X, Leonard EJ, et al. Phase II study of alisertib, a selective Aurora a kinase inhibitor, in relapsed and refractory aggressive B- and T-cell non-Hodgkin lymphomas. J Clin Oncol. 2014;32(1):44–50.
O'Leary B, Finn RS, Turner NC. Treating cancer with selective CDK4/6 inhibitors. Nat Rev Clin Oncol. 2016;13(7):417–30.
Rodon J, Brana I, Siu LL, De Jonge MJ, Homji N, Mills D, Di Tomaso E, Sarr C, Trandafir L, Massacesi C, et al. Phase I dose-escalation and -expansion study of buparlisib (BKM120), an oral pan-class I PI3K inhibitor, in patients with advanced solid tumors. Investig New Drugs. 2014;32(4):670–81.
Fritsch C, Huang A, Chatenay-Rivauday C, Schnell C, Reddy A, Liu M, Kauffmann A, Guthy D, Erdmann D, De Pover A, et al. Characterization of the novel and specific PI3Kalpha inhibitor NVP-BYL719 and development of the patient stratification strategy for clinical trials. Mol Cancer Ther. 2014;13(5):1117–29.
Campbell PM, Boufaied N, Fiordalisi JJ, Cox AD, Falardeau P, Der CJ, Gourdeau H. TLN-4601 suppresses growth and induces apoptosis of pancreatic carcinoma cells through inhibition of Ras-ERK MAPK signaling. J Mol Signal. 2010;5:18.
Wilhelm SM, Adnane L, Newell P, Villanueva A, Llovet JM, Lynch M. Preclinical overview of sorafenib, a multikinase inhibitor that targets both Raf and VEGF and PDGF receptor tyrosine kinase signaling. Mol Cancer Ther. 2008;7(10):3129–40.
Strumberg D, Schultheis B. Regorafenib for cancer. Expert Opin Investig Drugs. 2012;21(6):879–89.
Sierra RA, Trillo-Tinoco J, Mohamed E, Yu L, Achyut BR, Arbab A, Bradford JW, Osborne BA, Miele L, Rodriguez PC. Anti-jagged immunotherapy inhibits MDSCs and overcomes tumor-induced tolerance. Cancer Res. 2017;77(20):5628–38.
Meurette O, Mehlen P. Notch Signaling in the Tumor Microenvironment. Cancer Cell. 2018;34(4):536-48.
Shen Q, Cohen B, Zheng W, Rahbar R, Martin B, Murakami K, Lamorte S, Thompson P, Berman H, Zuniga-Pflucker JC, et al. Notch shapes the innate Immunophenotype in breast Cancer. Cancer Discov. 2017;7(11):1320–35.
Liu H, Wang J, Zhang M, Xuan Q, Wang Z, Lian X, Zhang Q. Jagged1 promotes aromatase inhibitor resistance by modulating tumor-associated macrophage differentiation in breast cancer patients. Breast Cancer Res Treat. 2017;166(1):95–107.
Fan X, Matsui W, Khaki L, Stearns D, Chun J, Li YM, Eberhart CG. Notch pathway inhibition depletes stem-like cells and blocks engraftment in embryonal brain tumors. Cancer Res. 2006;66(15):7445–52.
Manoranjan B, Venugopal C, McFarlane N, Doble BW, Dunn SE, Scheinemann K, Singh SK. Medulloblastoma stem cells: where development and cancer cross pathways. Pediatr Res. 2012;71(4 Pt 2):516–22.
Hoffman LM, Fouladi M, Olson J, Daryani VM, Stewart CF, Wetmore C, Kocak M, Onar-Thomas A, Wagner L, Gururangan S, et al. Phase I trial of weekly MK-0752 in children with refractory central nervous system malignancies: a pediatric brain tumor consortium study. Childs Nerv Syst. 2015;31(8):1283–9.
Kolb EA, Gorlick R, Keir ST, Maris JM, Lock R, Carol H, Kurmasheva RT, Reynolds CP, Kang MH, Wu J, et al. Initial testing (stage 1) by the pediatric preclinical testing program of RO4929097, a gamma-secretase inhibitor targeting notch signaling. Pediatr Blood Cancer. 2012;58(5):815–8.
Rayburn ER, Ezell SJ, Zhang R. Anti-inflammatory agents for Cancer therapy. Mol Cell Pharmacol. 2009;1(1):29–43.
Berkova Z, Tao RH, Samaniego F. Milatuzumab - a promising new immunotherapeutic agent. Expert Opin Investig Drugs. 2010;19(1):141–9.
Brown ER, Charles KA, Hoare SA, Rye RL, Jodrell DI, Aird RE, Vora R, Prabhakar U, Nakada M, Corringham RE, et al. A clinical study assessing the tolerability and biological effects of infliximab, a TNF-alpha inhibitor, in patients with advanced cancer. Ann Oncol. 2008;19(7):1340–6.
Pavlova NN, Thompson CB. The emerging hallmarks of Cancer metabolism. Cell Metab. 2016;23(1):27–47.
Beloribi-Djefaflia S, Vasseur S, Guillaumond F. Lipid metabolic reprogramming in cancer cells. Oncogenesis. 2016;5:e189.
Zhang J, Pavlova NN, Thompson CB. Cancer cell metabolism: the essential role of the nonessential amino acid, glutamine. EMBO J. 2017;36(10):1302–15.
DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv. 2016;2(5):e1600200.
De Santis MC, Porporato PE, Martini M, Morandi A. Signaling pathways regulating redox balance in Cancer metabolism. Front Oncol. 2018;8:126.
Martinez-Reyes I, Chandel NS. Mitochondrial one-carbon metabolism maintains redox balance during hypoxia. Cancer Discov. 2014;4(12):1371–3.
Kishton RJ, Rathmell JC. Novel therapeutic targets of tumor metabolism. Cancer J. 2015;21(2):62–9.
Akins NS, Nielson TC, Le HV. Inhibition of glycolysis and Glutaminolysis: an emerging drug discovery approach to combat Cancer. Curr Top Med Chem. 2018;18(6):494–504.
We thank all the contributors of Medulloblastoma advanced. Genomics International Consortium (MAGIC) project.
This study was supported by a grant from the National R&D Program for Cancer Control, Ministry for Health and Welfare, Republic of Korea (1420020), Mid-career Researcher Program through the National Research Foundation (NRF) grant funded by the Korea government (Ministry of Science and ICT, 2017R1A2B2008422), and NRF grants funded by the Ministry of Science, ICT & Future Planning, Republic of Korea (2012R1A1A1042953 and 2015R1A4A1041219). The funding bodies had no role in the design of the study, collection, analysis, and interpretation of data nor in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Goodness-of-fit tests for Cox proportional hazards model. Table S2. Subgroup-specific enriched gene sets (canonical pathways) and core genes positively and negatively associated with poor prognosis. Table S3. Cox regression analysis with chromosome aberrations occurred in more than 10% of patients in SHH subgroup. Table S4. Cox regression analysis with chromosome aberrations occurred in more than 10% of patients in Group 3 Table S5. Cox regression analysis with chromosome aberrations occurred in more than 10% of patients in Group4. Table S6. Cox regression analysis with chromosome aberrations in Group4 after excluding the patients with chr11 loss. Table S7. Cox regression analysis with chromosome aberrations in Group4 patients without chr11 loss nor i(17q). Table S8. Cox regression analysis with chromosome aberrations in Group4 patients harboring i(17q) but not chr11 loss. (XLSX 1264 kb)
Figure S1. Kaplan-Meier subgroup analysis. A, 10-year Kaplan-Meier subgroup analysis in the exploratory dataset. B, 10-year Kaplan-Meier subgroup analysis in the validation dataset. Figure S2. Association between the average relative expression of prognostic core genes and subtype and clinical characteristics. A, Boxplots of relative expression of 31 core genes from the top two PID pathways in SHH subgroup. B, Boxplots of relative expression of 20 core genes from the top two PID pathways in Group 3. C, Boxplots of relative expression of 31 core genes from the top two PID pathways in Group 4. A, large cell/anaplastic (LCA); C, Classic; D, Desmoplastic; M, Medulloblastoma with extensive nodularity (MBEN). *, Mann–Whitney U test p-value < 0.05; **, p-value < 0.01; ***, p-value < 0.001. Figure S3. Differential expressions of the prognostic genes identified in SHH subgroup. A, Boxplots of gene expressions in four subgroups. B, Differential gene expressions with Mann–Whitney U test p-values according to the status of chromosome aberrations in SHH subgroup. Figure S4. Differential expressions of the prognostic genes identified in Group 3. A, Boxplots of gene expressions in four subgroups. B, Differential gene expressions with Mann–Whitney U test p-values according to the status of chromosome aberrations in Group 3. Figure S5. Correlation between expressions of NOTCH and those of immune receptors. Scatter plots with regression lines are presented with Pearson’s correlation coefficients (r) calculated between the expressions of NOTCH1 (A), NOTCH2 (B), or NOTCH3 (C) and those of immune receptors that show r > 0.4. Figure S6. Differential expressions of the prognostic genes identified in Group 4. A, Boxplots of gene expressions in four subgroups. B, Differential gene expressions with Mann–Whitney U test p-values according to the status of chromosome aberrations in Group 4. (PDF 1629 kb)
About this article
Cite this article
Park, A.K., Lee, J.Y., Cheong, H. et al. Subgroup-specific prognostic signaling and metabolic pathways in pediatric medulloblastoma. BMC Cancer 19, 571 (2019). https://doi.org/10.1186/s12885-019-5742-x
- Signaling pathway
- Metabolic pathway
- Targeted therapy