3.1 Co-Expression Network Construction and Key Module Selection
To comprehensive understanding of antitumor immune responses in diffuse gliomas, we combined the cancer immunity-cycle score with diffuse gliomas RNA-seq profiles from the TCGA-LGG and TCGA-GBN database for a weighted gene co-expression network analysis. After removing outlying samples and filtering the genes with a MAD < 1, the expression profiles of 19,535 genes and 657 samples were extracted, and a sample clustering tree was drawn (Fig. 2A, B). When the soft threshold power value was set to 10, the scale independence reached 0.90, and as the scale-free topology criteria (Fig. 2C). When the minimum module size was set to 30 and the cut height was set to 0.25, all genes were divided into 24 different co-expression modules through the hybrid dynamic tree cutting, as shown in a cluster dendrogram (Fig. 2D). The members in each module were listed in Supplementary Material S2. Apart from the grey module consisting of many un-classified members, the yellowgreen module contained a minimum of 33 genes, while a maximum of 2,580 genes were included in the darkgrey module. Then, analyzing the correlation between each module and clinical phenotypes (Fig. 2E), which the darkgrey module had the highest positive correlation with the clinical phenotypes (Cor = 0.92, p<1e-200), while the greenyellow module had the highest negative correlation with the clinical phenotypes (Cor =0.8, p<1e-200) (Fig. 2F).
Fig. 2 Construction of WGCNA co-expression modules associated with the cancer-immunity cycle traits of diffuse gliomas and hub modules selection. A The cluster dendrogram of 669 patients with diffuse glioma samples. B Clustering dendrogram of 657 patients with diffuse gliomas and tumor immune cycle trait heatmap. C The scale-free topological model fit index [left] and the mean connectivity values (right) at different soft-thresholding powers. D Cluster dendrogram of genes based on a dissimilarity measure [1-TOM]. E Heatmap demonstrating the correlation between different module eigengene and the cancer immunity cycle. All genes were clustered into twenty-four modules, each color representing a type of module. Red represents a positive correlation, while blue represents a negative correlation. F The darkgrey module membership had the highest positive correlation with the cancer immunity cycle (Cor = 0.92, P < 1e-200), and the greenyellow membership module had the highest negative correlation with the cancer immunity cycle (Cor = 0.8, P < 1e-200).
Next, GO and KEGG pathway enrichment analyses of these two modules were performed through the DAVID database, respectively. Unsurprisingly, as shown in Fig 3A, the genes in the darkgrey module were significantly enriched in cell adhesion, inflammatory response, and angiogenesis of BPs; extracellular exosome, endoplasmic reticulum lumen, and extracellular region of CCs; protein binding, extracellular matrix structural constituent, and collagen binding of MFs; toxoplasmosis, focal adhesion, and ECM-receptor interaction of KEGG pathway. However, as shown in Figure 3B, genes in the greenyellow module were mainly enriched in chemical synaptic transmission, regulation of ion transmembrane transport, and neuron projection development of BPs; anchoring junction, glutamatergic synapse, and dendrite of CCs; ion channel binding, voltage-gated ion channel activity, and ionotropic glutamate receptor binding of MFs, glutamatergic synapse, synaptic vesicle cycle, and axon guidance of KEGG pathway.
Fig. 3 Functional enrichment analysis of genes in the darkgrey and greenyellow modules. A Left panel: GO enrichment analysis showed the top ten enriched BPs, CCs, and MFs in the darkgrey module. Right panel: KEGG enrichment analysis showed the top 30 enriched pathways in the darkgrey module. B Left panel: GO enrichment analysis showed the top ten enriched BPs, CCs, and MFs in the greenyellow module. Right panel: KEGG analysis showed the top 30 enriched pathways in the greenyellow module. BPs, biological processes; CCs, cellular components; MFs, molecular function
3.2 Identification of Hub Genes
Given that our study focused on immunosuppressive targets in diffuse gliomas, the greenyellow module, containing 1,128 genes was selected for further analysis. As described in the Methods, 95 hub genes (PPIhub) from the PPI network were identified from the PPI through the Degree algorithm of plug-in cytoHubba, with a connectivity degree ≥8. Additionally, 206 hub genes (MMhub) were identified in the greenyellow module based on high connectivity in the module connectivity threshold criterion |cor.geneModuleMembership| > 0.8. After taking the intersection of the Venn diagrams, the 37 common genes were identified among MMhub, and PPIhub (Fig. 4C). Subsequently, the PPI networks for the 37 common genes were constructed through the STRING database with a minimum interaction score ≥0.700. To figure out the potential biological processes and pathways of common genes, GO and KEGG pathway enrichment analyses were performed. As shown in Figure 4B, the BP was significantly enriched in the positive regulation of excitatory postsynaptic potential, chemical synaptic transmission, and exocytic insertion of neurotransmitter receptors to the postsynaptic membrane. The CCs were enriched in the glutamatergic synapse, anchoring junction, and Schaffer collateral - CA1 synapse. The significant MFs were mainly enriched in binding-related functions. In addition, KEGG Pathway analysis showed that they are associated with synaptic vesicle cycle, glutamatergic synapse, and cAMP signaling pathway (Fig. 4C). Through the degree algorithms of plug-in cytoHubba, the top 10 hub genes were calculated and were used for the final hub gene, including SNAP25, VAMP2, SYN1, DLG4, RIMS1, STXBP1, APBA1, CPLX1, PPFIA3, and UNC13A (Fig. 4D). Based on the GeneMANIA database, the co-expression network and related functions of hub genes were analyzed. These genes showed a complex PPI network with physical interactions of 48.64%, predicted of 20.82%, pathway of 14.98%, co-expression of 13.06%, and co-localization of 2.43% (Fig. 5E). GO analysis showed that hub genes are involved in the synaptic vesicle exocytosis, neurotransmitter secretion, and chemical synaptic transmission of BPs, the cytosol, glutamatergic synapse, and anchoring junction of CCs, and syntaxin-1 binding, syntaxin binding, and SNARE binding of MFs (Fig. 5F). These results emphasized the important role of synaptic-related functions in diffuse gliomas. In addition, KEGG Pathway analysis showed that six genes (RIMS1, SNAP25, UNC13A, STXBP1, CPLX1, and VAMP2) were involved in the synaptic vesicle cycle.
Fig. 4 Common gene expression profiles. A The MMhub and PPIhub showed an overlap of 37 common genes. B GO in which the common genes were involved. C KEGG pathway in which the common genes were involved. p-value < 0.05 was considered significant. D PPI network of the 37 common genes of diffuse gliomas. E Hub genes and their co-expression genes were analyzed using the GeneMANIA database. F GO enrichment analysis of the hub genes.
3.3 Validation of Hub Genes Expression
The progression of diffuse gliomas is usually accompanied by abnormal gene expression and poor prognosis. To verify the reliability of these gene expression levels. Microarray data sets containing the normal brain, LGG, and GBM were selected and the expression levels of these hub genes were analyzed. The results showed that compared with normal brain tissue, all hub genes [the expression value of APBA1 no statistically significant] were significantly down-regulated in LGG (Fig. 5A). Similarly, the expression of all genes in GBM was also lower than in normal brain tissues (Fig. 5B). Based on the GEPIA2 online tool, data from the match TCGA normal and GTEx tissue adjacent to the LGG and GBM tissues revealed a consistent trend (Fig. 6A). Moreover, survival analyses revealed that patients with diffuse gliomas in these genes had low expression groups as well as shorter OS and RFS outcomes in multiple datasets (Fig. 6B, C). Overall, our results suggested that these genes are lowly expressed in diffuse gliomas and are associated with a poorer prognosis.
Fig. 5 The expression level of the hub gene in GSE145372 dataset. The comparison between the two sets of data uses the mean T-test. P-value < 0.05 was considered statistically significant. LGG, low-grade gliomas; GBM, glioblastomas multiforme; Normal, control brain samples. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.
Fig. 6 The expression level and prognostic value of these genes in the GEPIA database. A The expression level of these genes comparison between the two sets of data based on the differential methods limma. |Log2FC|≥0.585 and q-value < 0.05 was considered to be significant. B, C Survival curves for prognostic analyses of these genes were evaluated using the OS and RFS method, respectively. The log-rank P<0.05 was applied to estimate statistical significance. LGG, low-grade gliomas; GBM, glioblastomas multiforme; N, Match TCGA normal and GTEx data; OS, overall survival; RFS, disease-free survival.
3.4 Immunological Correlation Analysis of Hub Genes
The immunological characteristics of hub genes are critical in determining diffuse gliomas that may benefit from antitumor immunotherapy. Our findings revealed that hub genes were negatively correlated with a large number of immunomodulators except for CCL27 deletion in diffuse gliomas (Fig. 6A). We also demonstrated that these hub genes were negatively correlated with the majority of immune checkpoints, including CD276, CD40, CD86, and IL2RB in diffuse gliomas (Fig. 6B) Next, the infiltration levels of TIICs in the TME were also assessed via the CIBERSORT algorithm. These hub genes were significantly negatively correlated with M2 macrophages, and positively associated with naive B cells, Plasma cells, Eosinophils, resting NK cells, follicular helper T cells, and activated NK cells (Fig. 6C). Furthermore, the hub genes in addition to were found to be negatively correlated with the activities of the majority of the steps in the cancer immunity cycle, including the release of cancer cell antigens [Step 1], and trafficking of immune cells to tumors (Step 4) (Th1 cell recruiting, Th2 cell recruiting, CD8 T cell recruiting, Basophil recruiting, Treg cell recruiting, CD4 T cell recruiting, Dendritic cell recruiting, Th22 cell recruiting, Macrophage recruiting, Monocyte recruiting, Neutrophil recruiting, NK cell recruiting, Eosinophil recruiting and MDSC recruiting), while were positively correlated with the cancer antigen presentation (Step 2), priming and activation (Step3), and killing of cancer cells (Step 7) (Fig. 6D). Taken together, these results suggest that the hub genes may be the potential unreported immunosuppressive target in diffuse gliomas.
Fig. 6 The effect of hub genes on immunological status in diffuse gliomas. A The correlation between the hub genes and 122 immunomodulators (the expression value of CCL27 missing in diffuse gliomas). B Correlation between the hub genes and 29 immune checkpoints. C The correlation between the hub genes and 28 tumor-associated immune cells was calculated with the CIBERSORT algorithm. D Correlation between the hub genes and steps of the cancer-immunity cycle. The color indicates the correlation coefficient. The asterisks indicate a statistically significant p-value calculated through Spearman correlation analysis. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.