Gene regulatory networks analysis for the discovery of prognostic genes in gliomas

Abstract

Gliomas are the most common and aggressive primary tumors of the central nervous system. Dysregulated transcription factors (TFs) and genes have been implicated in glioma progression, yet these tumors' overall structure of gene regulatory networks (GRNs) remains undefined. We analyzed transcriptional data from 989 primary gliomas in The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) to address this. GRNs were reconstructed using the RTN package which identifies regulons-sets of genes regulated by a common TF based on co-expression and mutual information. Regulon activity was evaluated through Gene Set Enrichment Analysis. Elastic net regularization and Cox regression identified 31 and 32 prognostic genes in the TCGA and CGGA datasets, respectively, with 11 genes overlapping, many of which are associated with neural development and synaptic processes. GAS2L3, HOXD13, and OTP demonstrated the strongest correlations with survival outcomes among these. Single-cell RNA-seq analysis of 201,986 cells revealed distinct expression patterns for these genes in glioma subpopulations, particularly oligoprogenitor cells. This study uncovers key GRNs and prognostic genes in gliomas, offering new insights into tumor biology and potential therapeutic targets.

Keywords: Glioma; Regulons; Survival.

Fig. 1

Least absolute shrinkage and selection operator (LASSO) analysis of regulons from TCGA and CGGA: (a) Schematic workflow. (b) Cross-validation results demonstrating the performance of LASSO for variable coefficient selection. The plot illustrates the partial likelihood deviance across different values of the regularization parameter (lambda), aiding in the selection of the optimal lambda. (c, d) LASSO coefficients for the identified regulons in the TCGA and CGGA dataset respectively. Each dot represents the magnitude and direction of the effect of each variable on the model. (e) Hierarchical tree-and-leaf representation of the identified regulons; the size of each circle corresponds to the number of genes within each regulon. Green circles indicate regulons derived from TCGA, while orange circles represent those from CGGA. (e) Heatmaps displaying the activity levels of the identified regulons across all samples from CGGA and TCGA.

Fig. 2

Characterization of genes within the regulons: (a) Functional enrichment analysis of all identified regulons from CGGA and TCGA, highlighting the biological processes significantly associated with the gene sets. (b) Venn diagram illustrating the overlap of genes within prognostically relevant regulons from CGGA and TCGA. The diagram indicates the number of unique and shared genes. (c, d) Cross-validation results for Elastic Net variable coefficient screening, demonstrating the stability and performance of selected variables across different parameter settings. The plots shows partial likelihood deviance as a function of the regularization parameters, aiding in the identification of the optimal model, significant variables coefficient values and elastic net coefficients for the identified regulons in the CGGA and TCGA dataset respectively. Each dot represents the magnitude and direction of the effect of each variable on the model.

Fig. 3

Common prognostic genes. (a) Upset plot illustrating the overlap of genes identified by elastic net in the TCGA and CGGA datasets. (b) Hazard ratios for the eleven common prognostic genes, adjusted for tumor grade and patient age. The plot presents the hazard ratios with corresponding confidence intervals. (c, d) Survival curves of CGGA and TCGA respectively, stratified by the maximization of the log-rank test statistic, demonstrating the differential survival probabilities of patients with high versus low expression levels. The curves are accompanied by log-rank test p-values to assess the statistical significance of the observed differences. Benjamini–Hochberg (BH) was used to adjust for multiple hypothesis.

Fig. 4

Gene expression at the single-cell level. (a) Schematic representation of the analysis workflow used to assess gene expression at the single-cell level. (b) UMAP plot displaying the distribution of all 201,989 single cells, with cells color-coded according to their respective labels. (c) Feature plots illustrating the expression levels of the eleven prognostic genes across all cells. Each plot showcases the spatial distribution of gene expression within the cell population. (d) Dotplot showcasing average expression (color gradient) and percentage of cells expressing (dot size) of genes across cell populations.

Fig. 5

Glioma subcluster exploration. (a) UMAP plot illustrating glioma subclusters identified through clustering analysis. Seven distinct GCs are specified. (b) UMAP plot depicting the cell cycle phases of the cells included in the analysis. Cells are color-coded according to their respective cell cycle phases (G1, S, G2/M). (c) Heatmap displaying the top 5 marker genes for each GC. (d) UMAP plot showcasing the module scores of the four glioma cell states (AC, MES, NPC, OPC). (e) Quadrant plot illustrating the distribution of each GCs across four quadrants, each representing one of the glioma cell states.

Fig. 6

Gene expression across glioma subclusters (GCs). (a) Feature plot displaying the expression levels of eleven selected genes within the glioma subset. Each cell is color-coded according to the expression intensity of the respective gene. (b) Dotplot of the mean and percent expression from the eleven genes across GCs.

Fig. 7

Proposed mechanisms linking prognostic genes to glioma biology: Schematic processes that each prognostic gene is involved in the glioma biology.