3.1 Difference of immune cell infiltration among medulloblastoma subgroups.
Tumor-infiltrating immune cells in medulloblastoma samples were examined with MethylCIBERSORT package in the R software based on the methylation data [12]. Based on the transcriptome expression data of tumor samples, estimate score, immune score and stromal score are estimated by using estimate according to the proportion of stromal and immune cells, which are used to predict tumor purity. The higher the estimated score, the lower the purity of tumor. The degree of immune cell infiltration, the immune scores, and purity of tumors of 763 samples in GSE85212 are displayed in Fig. 1A. It was evaluated by transcriptomic data within the same group of the samples from GSE85212. Variance analysis was employed to compare the degree of immune cell infiltration among different subtypes of MB and to establish the differential infiltrating immune cells. Results indicated that among the four subgroups of MB, 9 kinds of infiltrating immune cell (such as CD14, Treg, CD58, CD8, CD19, etc) were statistically significant in the degree of tumor infiltration (Figure. 1B) (p < 0.001).
3.2 Differential immune expressed gene and functional analysis among medulloblastoma subgroups.
Merging the immune gene sets from ImmPort and InnateDB databases, in total 2533 immune-related genes, consisting of 1794 immune gene from ImmPort immune gene set and 1052 immune gene from InnateDB immune gene set were included. Limma, a package for analyzing the differential expression, was utilized to screen the differentially expressed immune genes (imm-DEGs) among four subgroups of medulloblastoma. 293 immune genes that expressed differentially among four subgroups of medulloblastoma was established (adjusted p Value < 0.05). The differential gene distribution and heatmap of differential expression genes are illustrated in Fig. 2A-B.
The biological function enrichment analysis, including GO analysis and KEGG pathway enrichment analysis, was utilized. The GO analysis results emphasized that 293 imm-DEGs were associated with regulation of chemotaxis, regulation protein serine/threonine kinase activity, MAP kinase activity in the BP category(Figure.2C); receptor ligand activity, growth factor binding and cytokine receptor binding in the MF category(Figure.2D); and transcription factor complex, nuclear transcription factor complex, membrane microdomain in the CC category(Figure.2E); in the KEGG pathways, including RAS signaling pathway, Human cytomegalovirus infection, Chemokine signaling pathway, T cell receptor signaling pathway, etc. were all implicated (Figure.2F).
3.3 Construction of the co-expression network based on imm-DEGs and identify key modules.
WGCNA, based on the imm-DEGs expression profiles, was employed to construct the co-expression network and distinguish the gene with similar expression pattern. A hierarchical clustering (Figure.3A) was used by setting the soft threshold power to 3, cut height to 0.25, and the minimum module size to 30, to cluster the gene into different modules with different reprehensive color. In the end, four gene modules were established (Figure.3B).
The degree of immune cell infiltration of samples was collected as a trait, and module-trait relationships were calculated according to the correlation in order to match the modules to their strongly related traits. As shown in Fig. 3C, the blue module was significantly related to Fibroblast cells, Eos cells and Endothelial cells, which are important components of immune microenvironment, and their interaction with tumor cells plays an important role in the growth of cancer (corr = 0.65, 0.55, -0.6, p value = 5e-92, 1e-60, 3e-76; respectively) and the correlation was stronger than other modules. Thus, blue was chosen to be the key module.
3.4 Hub genes recognition.
Genes from the selected blue key module were used to construct the PPI network, in which 18 hub genes were established by setting the node degree value ≥ 5(Figure.4A). 25 candidate hub genes from the blue module were identified by setting the gene significance (GS) > 0.4 and module membership (MM) > 0.6 (Figure.4B). Three immune hub genes (GAB1, ABL1, CXCR4) were obtained by synthesizing the hub genes in PPI network and the candidate hub genes in blue module. (Figure. 4C).
The corrplot package in the R software was employed to investigate the correlation within the 3 hub genes. Results indicated that hub gene GAB1 was positively correlated with CXCR4 (Pearson corr = 0.77); ABL1 was negatively correlated with both CXCR4 and GAB1 (Pearson corr =-0.55, -0.47, respectively); these three hub genes were related with each other (Figure.4D).
The correlation between hub genes and the differential immune infiltrating cells were also scrutinized. As the result shows, both hub genes GAB1 and CXCR4 have a significant positive correlation with CD4_Eff and Fibroblast cells (p < 0.01), and a negative correlation with Endothelial cells and the Eos cells (p < 0.001); while hub gene ABL1 has an opposite trend (Figure.4E). In addition, we selected the characteristic markers of these correlation immune cells, the immune hub genes and immune cells were confirmed by immunohistochemistry (Fig. 6B-C and Supplementary Figure.1). Due to the limited number of tissue samples, we will accumulate more samples for verification in the future.
3.5 Medulloblastoma subgroups markers recognition based on hub genes.
Three hub genes expression patterns in the subgroups of medulloblastoma were also studied and markers of the subgroups were recognized. It was indicated that immune hub gene GAB1 and CXCR4 express higher in the SHH subgroup than other three subgroups; immune hub gene ABL1 was expressed weakly in the SHH subgroup, while expressed strongly in other three subgroups (Figure.5A). The methylation level of immune hub genes among four subgroups was established by combining the immune hub gene with the DNA methylated data. The results of one-way ANOVA showed that the methylation levels of GAB1 and CXCR4 were significantly different among the four subtypes (p < 0.001) (Figure.5B).
3.6 Further verification of hub genes in independent data set and at histological level.
GSE37418 was downloaded from GEO databases, which contained 73 medulloblastoma samples. GSE37418 consisted of 10 samples of SHH subgroup, 8 samples of WNT subgroup, 16 samples of Group 3 subgroup, and 39 samples of Group 4. Comparison and analysis of the expression patterns of subgroup markers were achieved. We found that the expression of GAB1 and CXCR4 was higher in SHH subtype than the other three subtypes, while ABL1 expressed lower in SHH subtype and higher in the other three subtypes (Figure.6A). Next, the immune hub genes were confirmed by immunohistochemistry (Fig. 6B-C). The expression of CXCR4 and GAB1 was apparently high in SHH subgroups compared with other subgroups. As Fig. 6B illustrates, the expression of ABL1 was the lower compared with CXCR4 and GAB1. Results of immunochemistry staining were consistent with the result of our bioinformatic analysis.
3.7 A multi-factor regulatory network of medulloblastoma based on immune hub genes.
A multi-factor regulatory network consisted of miRNA, lncRNA, and transcriptional factor (TF) was designed, combined with public databases (such as starBase, Harmonizome) and immune hub genes. Interactions among genes and their following products play a crucial role in many biological processes [13]. In the present study, a multi-factor network regarding hub genes combined with public databases (such as starBase, Haemonizome) was designed and constructed, which included miRNA, lncRNA, mRNA, and TF. Figure 7A shows the multifactor regulatory network of immune hub gene ABL1, GAB1 and CXCR4 in medulloblastoma. We found that CXCR4 is mainly involved in the regulation of lncRNA. GAB1 is related to the regulation of miRNA. ABL1 is not only involved in the regulation of lncRNA and miRNA, but also in the regulation of transcription factors. MiRNA, lncRNA and TF participate in the regulation of gene expression at the transcriptional and post transcriptional levels, and play a crucial role in gene expression. Therefore, they can be used as potential drug targets to indirectly regulate the expression of hub gene by acting on miRNA, lncRNA or TF.