Identification of DEGs
Following differential gene expression analysis and filtering with the cut-off criteria of P. value <0.05 & |logFC|>1, the DEGs between tumor and normal cells were screened. For the GSE43290 dataset, a total of 1499 DEGs were identified, out of which 73 were up-regulated and 1426 were down-regulated. For the GSE77259 dataset, 4830 DEGs, including 2536 up-regulated and 2294 down-regulated genes, were identified. Moreover, a total of 673 overlapping DEGs, such as MYC, TNFAIP3 and SLC2A3 in both GSE43290 dataset and validation GSE77259 dataset were obtained (Supplementary figure 1).
Furthermore, volcano plot and PCA clustering of DEGs in GSE43290 showed that the differentially expressed DEGs between the two groups could be significantly distinguished (Figure 1A and 1B).
Furthermore, using a threshold of |logFC|>2, 255 DEGs were analyzed using the bilayer cluster analysis (Figure 1C) to further predict the potential functions of the dysregulated mRNAs in meningioma, detected by enrichment analysis. The results of GO analysis showed that these enriched mRNAs were functionally classified by Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) (Figure 1D). Among them, several BPs seemed to be related to the mechanisms of HF; moreover, CCs and MFs were found to be mainly associated with regulation of cardiac electrical activity. Furthermore, KEGG pathway analysis revealed a series of enriched pathways, such as IL-17 signaling pathway and TNF signaling pathway (Figure 1E)
The immune cell infiltration landscape
Based on the gene expression profiles of GSE43290, the difference in immune infiltration between tumor and control samples among 22 immune cells types, were detected using CIBERSORT algorithm. As shown in Figure 2A, 48 of the 51 samples were valid (P value <0.05), including 44 tumor samples and 4 normal samples. Among the 22 immune cell types, the percentage of CD8+ T cells (yellow) and M2 macrophages (blue) in each sample was relatively high. The heatmap showing the 22 tumor immune cells was illustrated in Figure 2B. Based on the violin plot (Figure 2C), the probability of distribution of memory B cells, regulatory T cells (Tregs) and resting mast cells in tumor were obviously higher than those in the control samples (all, P<0.05), indicating that the immune infiltration proportions could be useful in distinguishing meningioma patients from healthy patients.
In addition, a lot of literature have shown that mast cells (MC) were associated with the development of meningioma [13, 14, 27], and resting mass cells infiltration ratio was relatively high in this study. Hence, the follow-up analysis focused on resting mast cells.
Weighted correlation network analysis (WGCNA) co-expression networks
When correlation coefficient threshold was set at 0.85, the soft-thresholding power was 6 (Supplementary figure 2A). Through WGCNA analysis, 10 co-expression modules were constructed based on the 1499 DEGs (Supplementary figure 2B), of which the grey module was the gene set that could not gather into other modules, so there were only 9 valid gene modules.
Subsequently, as shown in Figure 3C, the correlation between each module and the degree of invasion of immune cells in tumor tissue was calculated. The results revealed that yellow, the module with the highest correlation, contained 158 DEGs, consisting of 2 up-regulated and 156 down-regulated genes (Supplementary table 1). Yellow modules are negatively related to degree of resting mast cells infiltration in meningioma tissues, and the correlation between genes in the yellow module and the degree of immune cell infiltration was summarized in the Supplementary table 2.
Furthermore, GO functional analysis revealed that the genes in yellow module were dominant in biological processes involved in response to leukocyte migration, molecules of bacterial origin, and inflammatory response (Figure 4A); KEGG pathway analysis revealed that the genes were mainly enriched in the TNF signaling pathway, cytokine-cytokine receptor interaction and IL-17 signaling pathway (Figure 4B). The path diagram was annotated based on the pathview package in R, (Figure 4C).
Risk model construction of key genes
In order to identify whether each key gene related to mast cells resting correlated with prognosis of meningioma patients, univariate and multivariate analyses were performed. K-M survival analysis revealed that CXCL8 and MYC were associated with prognosis of meningioma patients. Moreover, based on the COX univariate and multivariate regression analyses showed CXCL8 and MYC might be prognostic biomarkers for meningioma patients (Table 1 and 5A).
Subsequently, based on the 9 genes in univariate regression analysis above, the risk model was constructed. Firstly, the meningioma patients were categorized into low- or high-risk patients, and the cut off value was the median risk score (Figure 5B). The scatter plot of survival time of risk model samples showed that the survival time of samples was relatively lower in the high-risk group the low-risk group (Figure 5C). Monolayer clustering analysis of RNA expression in the risk model of each sample showed differences in the key genes between high and low risk groups (Figure 5D). Moreover, K-M survival curve indicated that the survival time between the high and low risk groups was significantly different (P value = 0.00436; Figure 5E), suggesting the risk model of 9 key gene related to mast cells resting was successfully constructed.
TF-miRNA- mRNA co-regulation network
Based on the database (miRWalk3.0, TargetScan, MiRDB, and MirTarBase) with the Score > 0.95, 145 miRNAs acting on 3'UTR region of the genes in the risk model associated with mast cells-resting immune cells were predicted. Then, according to the retrieval of HMDD V3.2 database, 3 miRNA-mRNA pairs (including MYC-miR-145-5p, TNFAIP3-miR-29c-3p and TNFAIP3-miR-335-3p) were obtained, respectively. Subsequently, and 37 TF-mRNA pairs (including 6 mRNAs and 27 TFs) were screened using the online database of TRRUST. Lastly, based on mRNA-miRNA pairs and TF-mRNA pairs, miRNA-TF-mRNA co-regulation network (ceRNA) was constructed using Cytoscape (Figure 6A).
Additionally, in order to further investigate the potential biological functions of the genes involved in ceRNA network, the GO terms and KEGG pathway analyses were conducted. As illustrated in Figure 6B, the functional processes of the genes were primarily related to biological regulation, while KEGG pathway analysis revealed that the genes were enriched in Tumor Necrosis Factor (TNF) signaling pathway and IL-17 signaling pathway (Figure 6C).
PPI networkand modules analysis
According to STRING database, the PPI network (PPI score = 0.4) was constructed with 27 TFs in miRNA-TF-mRNA co-regulation network and 9 DEGs in the risk model using Cytoscape. As shown in Figure 7A, there were 36 nodes and 197 interactions in the PPI network. With score > 12, a module with 15 nodes and 88 interactions was further revealed from the PPI network using the MCODE of Cytoscape software (Figure 7B).
Drug-gene network
Similarly, based on the 9 DEGs in the risk model, a total of 50 drug-gene interacting pairs, including 3 target genes(CXCL2, CXCL8 and MYC), and 49 drugs were identified using the DGIdb 3.0 database. Furthermore, drug-gene interaction network was constructed by Cytoscape software (Figure 8).