Identification of microglia-associated hub genes and a CX3CR1/P2RY12 two-gene diagnostic model for cross-subtype epilepsy-associated cortical malformations: integrated transcriptomic analysis with external validation

1 Introduction

Epilepsy affects approximately 50 million people worldwide, with approximately 30% of patients developing resistance to anti-seizure medications despite optimal pharmacotherapy (Thijs et al., 2019). Focal cortical dysplasia (FCD) and tuberous sclerosis complex (TSC) represent the most common structural etiologies encountered in epilepsy surgery series (Blumcke et al., 2017; Najm et al., 2022). FCD is classified histopathologically into types I, IIa, IIb, and III, characterized by cortical dyslamination, dysmorphic neurons, and balloon cells (Blumcke et al., 2017). TSC, caused by loss-of-function mutations in TSC1 or TSC2, produces cortical tubers that share histological features with FCD type IIb (Curatolo et al., 2015). Despite distinct genetic etiologies, the clinical and histopathological overlap suggests shared downstream molecular mechanisms.

Neuroinflammation – particularly microglial activation – is increasingly recognized as a contributor to epilepsy pathogenesis (Eyo et al., 2017; Vezzani et al., 2019). Activated microglia have been demonstrated in surgical specimens from FCD and TSC patients, and TSPO-PET imaging has revealed microglial activation in epileptogenic foci (Gershen et al., 2015). Pro-inflammatory cytokines including IL-1beta, TNF-alpha, and IL-6 are elevated in cerebrospinal fluid and brain tissue from epilepsy patients (Vezzani et al., 2019). However, the transcriptome-wide molecular programs driving microglial activation across different epilepsy subtypes remain incompletely characterized.

Systems biology approaches, including weighted gene co-expression network analysis (WGCNA) and machine learning-based feature selection, offer powerful frameworks for identifying disease-associated gene modules and diagnostic biomarkers from high-dimensional transcriptomic data (Friedman et al., 2010; Langfelder and Horvath, 2008). Gene set enrichment analysis (GSEA) provides complementary information to traditional over-representation analysis by considering the rank order of all genes, enabling detection of coordinated pathway-level changes (Wu et al., 2021). Immune cell deconvolution methods such as MCP-counter and xCell enable characterization of the tissue immune microenvironment from bulk RNA-seq data (Aran et al., 2017; Becht et al., 2016).

In this study, we primarily focused on FCD and TSC as the most common cortical malformations, while also including TLE-HS as a comparator epileptogenic lesion to identify broader shared molecular signatures. Specific aims were: (1) identify shared DEGs across all four subtypes; (2) construct WGCNA co-expression networks and identify cross-subtype hub genes; (3) build and externally validate a multi-gene diagnostic classifier using LASSO and random forest; (4) characterize immune microenvironment alterations using multi-method deconvolution; and (5) perform GSEA, drug repurposing, and transcription factor prediction. In addition, we performed in vitro validation experiments to confirm the protein-level expression patterns of key hub genes under disease-relevant stimulatory conditions.

2 Materials and methods

2.1 Data acquisition and preprocessing

Transcriptomic data were obtained from the Gene Expression Omnibus (GEO)¹. The training dataset GSE256068 comprised RNA-seq data from 85 samples in total, including 64 cortical samples (FCD2a, n = 18; FCD2b, n = 18; TSC, n = 14; ControlCortex, n = 14) and 21 hippocampal samples (TLE-HS, n = 8; ControlHippocampus, n = 13), profiled on the Illumina HiSeq platform. Clinical metadata included tissue region, sex, age, and seizure frequency (Supplementary Table 1). The external validation dataset GSE213488 included RNA-seq count data from cortical gray matter samples (FCD2a, n = 7; FCD2b, n = 6; Control, n = 4). Raw count matrices were normalized to log2(CPM+1) values prior to analysis. No ethical approval was required as all data were obtained from publicly available repositories. Principal component analysis of cortical samples is shown in Supplementary Figure 7.

2.2 Differential expression analysis

Differential expression analysis was performed using the limma package (v3.60) in R (v4.5.1) (Ritchie et al., 2015). For cortical samples, the linear model included group (ControlCortex, FCD2a, FCD2b, TSC) and tissue region (temporal, frontal, occipital, parietal) as covariates. For hippocampal samples, the model was ∼ group. Contrasts were constructed for FCD2a vs. ControlCortex, FCD2b vs. ControlCortex, TSC vs. ControlCortex, and TLE-HS vs. ControlHippocampus. DEGs were defined by |log2 fold change| > 1 and Benjamini-Hochberg adjusted P-value < 0.05. Gene IDs were mapped to HGNC symbols using org.Hs.eg.db (v3.19). Intersection analysis was performed using UpSetR and ggVennDiagram (Gao et al., 2021), and volcano plots were generated. As a sensitivity analysis, a three-way intersection restricted to cortical samples (FCD2a, FCD2b, and TSC), excluding hippocampal TLE-HS, was additionally performed (Supplementary Figure 6), which identified 564 shared DEGs across the three cortical groups. All 49 four-way shared DEGs remained shared among the three cortical groups, and all 15 hub genes and all 15 KEGG pathways were retained, confirming that the core immune-related signature is robust to the exclusion of TLE-HS.

2.3 Weighted gene co-expression network analysis

Weighted gene co-expression network analysis was performed using the WGCNA package (v1.72) (Langfelder and Horvath, 2008). The top 5,000 most variable genes (by median absolute deviation) from cortical samples were used to construct a signed co-expression network with soft-thresholding power beta = 4 (scale-free R-squared > 0.85). Modules were identified using dynamic tree cutting (deepSplit = 2, minModuleSize = 30) and merged at a cut height of 0.25. Module eigengenes were correlated with clinical traits. The complete module gene assignments are provided in Supplementary Table 3. Module membership (MM) and gene significance for TSC (GS.TSC) were calculated using TSC as a representative disease phenotype because TSC is driven by a well-defined monogenic pathway (TSC1/TSC2-mTORC1) and showed the strongest association with the brown module. Importantly, this choice does not bias the final gene set toward TSC-specific genes: hub gene candidates were additionally required to overlap with the 49 DEGs shared across all four subtypes, ensuring cross-subtype differential expression. Gene significance for FCD2a and FCD2b was also computed for the 15 hub genes (Supplementary Table 4). Hub genes were defined by: (1) |GS.TSC| > 0.5; (2) |MM| > 0.8 in brown or red modules; (3) membership in the 49 shared DEGs. In silico node deletion analysis (also termed virtual knockout) was performed by computing mean connectivity of the brown module TOM before and after removing each gene.

2.4 Functional enrichment analysis

Over-representation analysis (ORA) of Gene Ontology (GO) and KEGG pathways was performed using clusterProfiler (v4.16) with Benjamini-Hochberg correction (adjusted P < 0.05) (Wu et al., 2021). Gene set enrichment analysis (GSEA) was performed using clusterProfiler::gseKEGG and gseGO with gene lists ranked by log2 fold change for each comparison. KEGG pathways significantly enriched (adjusted P < 0.05) in all four comparisons were identified as cross-subtype common pathways (Supplementary Tables 5, 6).

2.5 Protein-protein interaction network

A protein-protein interaction (PPI) network was constructed using the STRING database (v12.0) with a minimum interaction score of 0.4 (Szklarczyk et al., 2023).

2.6 Immune cell infiltration analysis

Immune cell infiltration was assessed using MCP-counter (Becht et al., 2016) and independently cross-validated using xCell (Aran et al., 2017). Differences across groups were assessed by Kruskal-Wallis tests, and associations with clinical variables were evaluated by Spearman correlation and partial correlation analysis (ppcor package).

2.7 Diagnostic model construction and validation

The 15 identified hub genes were used as the input feature set for machine learning-based feature selection. LASSO logistic regression (glmnet v4.1, lambda.min, 5-fold CV) (Friedman et al., 2010) and random forest (randomForest v4.7, 500 trees) (Breiman, 2001) were applied to identify the most discriminative features. For random forest, variable importance was assessed by mean decrease in Gini impurity with ntree = 500 and default mtry (square root of the number of features). Genes selected by both algorithms were used to construct a multi-gene logistic regression model. Performance was assessed by ROC-AUC (pROC v1.19) (Robin et al., 2011). To avoid information leakage, feature selection was performed within each cross-validation fold using a nested approach: within each training fold, LASSO and random forest were applied to select features, and the resulting model was evaluated on the held-out fold. Across the five folds, CX3CR1 was selected in all five folds and P2RY12 in four of five folds, supporting their robustness as diagnostic features. The model was subsequently evaluated by internal 5-fold cross-validation and external validation on GSE213488. For cross-platform validation, expression values in both datasets were log2(CPM+1)-transformed to ensure comparability. No explicit batch correction was applied, and this is acknowledged as a limitation.

2.8 Drug repurposing and transcription factor prediction

Drug repurposing and transcription factor enrichment were performed using Enrichr (Chen et al., 2013) with the 15 hub genes as input, querying the DSigDB, DrugMatrix, ChEA_2022, and TRRUST databases. miRNA target prediction was performed using TargetScan and miRTarBase. Terms with adjusted P < 0.05 were considered significant.

2.9 Statistical analysis

All analyses were performed in R v4.5.1 with Bioconductor 3.20. Multiple testing was corrected using the Benjamini-Hochberg method. Group comparisons used Kruskal-Wallis or Wilcoxon rank-sum tests. Correlations were assessed using Spearman’s rho. Partial correlation controlled for age, sex, and disease subtype. P < 0.05 was considered significant.

For in vitro experiments, comparisons among the three groups (Control, LPS, MHY1485) were performed using one-way ANOVA followed by Tukey post-hoc test. Data are presented as mean ± SEM from three independent experiments. P < 0.05 was considered statistically significant.

2.10 Cell culture and in vitro validation

BV-2 murine microglial cells were cultured in DMEM supplemented with 10% fetal bovine serum at 37 °C with 5% CO₂. Cells were divided into three groups (n = 3 per group): (1) Control group, treated with phosphate-buffered saline (PBS); (2) LPS group, stimulated with lipopolysaccharide (100 ng/mL, 24 h) to induce classical neuroinflammation; and (3) MHY1485 group, treated with MHY1485 (10 μM, 24 h), a potent mTORC1 activator (Laplante and Sabatini, 2012), to mimic the mTOR pathway hyperactivation characteristic of FCD and TSC (Weichhart et al., 2008).

TABLE 1. Top 20 differentially expressed genes per comparison.

Total RNA was extracted using TRIzol reagent (Invitrogen) and reverse-transcribed into cDNA using the PrimeScript RT reagent kit (TaKaRa). Quantitative real-time PCR (qPCR) was performed using SYBR Green Master Mix (TaKaRa) on an Applied Biosystems QuantStudio 6 Flex system. Primer sequences, sourced from PrimerBank and published literature, were as follows: Irf8 (Forward: 5′-CAGGTGGGCAAGATCCAGAA-3′, Reverse: 5′-TGGTGAGGAAGTTGGCGATA-3′), Cx3cr1 (Forward: 5′-CCATCTGCTCAGGACCTCAC-3′, Reverse: 5′-CA CCAGACACTCGCCATTTCT-3′), P2ry12 (Forward: 5′-TGCCAT CTGTACCTGCCTTC-3′, Reverse: 5′-CGTCCCTTTTCCCCAT CCAA-3′), Tlr7 (Forward: 5′-GGAAATTGCCCTCGATGTTA-3′, Reverse: 5′-CAAAAATTTGGCCTCCTCAA-3′), and Gapdh (Forward: 5′-AGGTCGGTGTGAACGGATTTG-3′, Reverse: 5′-TGTAGACCATGTAGTTGAGGTCA-3′). Relative mRNA expression was calculated using the 2 to the power of negative delta Ct method with Gapdh as the internal reference.

For Western blot analysis, total protein was extracted using RIPA lysis buffer (Beyotime) containing protease and phosphatase inhibitors, separated by 10% SDS-PAGE, and transferred to PVDF membranes (Millipore). Membranes were blocked with 5% non-fat milk and incubated overnight at 4 °C with the following primary antibodies: anti-IRF8 (1:1000, Proteintech, 18977-1-AP), anti-CX3CR1 (1:500, Proteintech, 29819-1-AP), anti-P2RY12 (1:500, Proteintech, 11976-1-AP), anti-phospho-S6 (Ser240/244) (1:1000, Proteintech, 29220-1-AP), and anti-GAPDH (1:5000, Proteintech, 60004-1-Ig). After washing, membranes were incubated with HRP-conjugated secondary antibodies (1:5000, Proteintech) for 1 h at room temperature. Protein bands were visualized by enhanced chemiluminescence (ECL, Beyotime) and quantified using ImageJ. All experiments were performed using three independent biological replicates (n = 3), each with technical duplicates.

3 Results

3.1 Forty-nine shared DEGs across four epilepsy subtypes

Differential expression analysis identified substantial transcriptomic alterations: 865 DEGs (FCD2a vs. Control), 1,915 DEGs (FCD2b vs. Control), 2,772 DEGs (TSC vs. Control), and 3,749 DEGs (TLE-HS vs. Control) (Figures 1A–D, Supplementary Figure 1, and Table 1). The majority of DEGs were upregulated. Four-way intersection analysis revealed 49 shared DEGs across all subtypes (Figures 1E, F and Supplementary Table 2), all consistently upregulated, indicating a conserved activation signature.

3.2 WGCNA reveals disease-associated immune modules

Weighted gene co-expression network analysis of cortical samples (beta = 4, signed network) identified the brown module (MEbrown) and red module (MEred) as the modules most significantly associated with disease status (Figure 2 and Supplementary Figure 3). Gene Ontology enrichment of brown module genes revealed predominant involvement in immune system processes, inflammatory response, and myeloid leukocyte activation.

3.3 Fifteen hub genes converge on microglial identity

Hub gene screening by GS, MM, and DEG overlap identified 15 hub genes, all functionally associated with microglia/myeloid lineage: IRF8, CLEC7A, CYBB, TLR7, CX3CR1, CCR1, C3AR1, P2RY12, IRF5, CD300C, GPR34, MPEG1, RGS18, LILRA2, and CD84 (Figure 3A–C; Table 2). CX3CR1 and P2RY12 are canonical homeostatic microglial markers (Hickman et al., 2013). IRF8 is a master transcription factor for microglial development (Kierdorf et al., 2013). TLR7 functions as an endosomal innate immune sensor. C3AR1 mediates complement C3a signaling (Stevens et al., 2007). IRF5 and CYBB drive inflammatory effector functions (Takaoka et al., 2005). This highly non-random functional convergence implicates a conserved microglia/myeloid-associated program in epileptic cortex.

TABLE 2. Fifteen hub genes identified by WGCNA and three-way DEG overlap.

3.4 GSEA reveals convergent innate immune pathway enrichment

GSEA identified 15 KEGG pathways significantly enriched across all four comparisons (Figure 4, Table 3, and Supplementary Figures 2, 4, 5). The shared pathways clustered into six biological themes: (1) innate immune sensing (Toll-like receptor, NOD-like receptor, C-type lectin receptor signaling); (2) inflammatory signal transduction (NF-kB, JAK-STAT, TNF signaling); (3) immune cell communication (cytokine-cytokine receptor interaction, chemokine signaling); (4) phagocytosis (phagosome, Fc-gamma-R-mediated phagocytosis); (5) complement and coagulation cascades; and (6) hematopoietic/immune cell biology (hematopoietic cell lineage, osteoclast differentiation, Th17 differentiation). Notably, three hub genes directly participate in these pathways: TLR7 in TLR signaling, IRF5/IRF8 in NF-kB signaling, and CCR1/CX3CR1 in chemokine signaling.

TABLE 3. Fifteen common KEGG pathways significantly enriched across all four comparisons (NES and adjusted P-value).

NES, normalized enrichment score from gene set enrichment analysis. *P < 0.05, **P < 0.01, ***P < 0.001 (adjusted P-value, Benjamini-Hochberg method). Negative NES values indicate downregulation in the disease group.

TABLE 4. Per-gene AUC values across training comparisons and external validation.

3.5 Myeloid and T-cell infiltration correlates with clinical features

MCP-counter analysis revealed significantly elevated Monocytic lineage and T cell scores in FCD2a, FCD2b, and TSC compared to ControlCortex (Figure 5A). xCell independently confirmed monocyte/macrophage and dendritic cell expansion (Figure 5D). Full xCell immune scores are provided in Supplementary Table 7. Exploratory analyses showed that monocytic lineage scores showed a non-significant positive trend with seizure frequency in FCD2a (Spearman rho = 0.478, P = 0.052; Figure 5B). No correction for multiple comparisons was applied to these exploratory immune-cell/clinical-variable correlation analyses, and T cell scores showed a significant negative correlation with patient age (Spearman rho = −0.373, P = 1.34e-3; Figure 5C).

3.6 PPI network and virtual knockout analysis

STRING analysis revealed a dense PPI network centered on an IRF8-TLR7-CX3CR1 axis (Figure 6A). While the PPI network highlights the central connectivity of the IRF8-TLR7-CX3CR1 axis, in silico node deletion analysis additionally identified CLEC7A, LILRA2, and IRF5 (followed by C3AR1 and CD84) as the hub genes whose removal caused the largest decreases in brown module connectivity (Figure 6B), providing independent evidence of their structural importance.

3.7 CX3CR1 + P2RY12 two-gene diagnostic model

Individual per-gene ROC curves for all 15 hub genes are shown in Figure 7A. LASSO selected five genes (CX3CR1, C3AR1, P2RY12, IRF5, RGS18), and random forest ranked CX3CR1, P2RY12, GPR34, CD84, and CD300C as the top five features (Figures 7B, C). Algorithm intersection yielded CX3CR1 and P2RY12. A logistic regression model incorporating these two genes achieved an AUC of 0.979 (95% CI: 0.947–1.000) in the training set (Figure 7D). Nested five-fold cross-validation, with feature selection performed within each training fold, yielded an AUC of 0.862 (95% CI: 0.693–1.000), an accuracy of 94.12%, a sensitivity of 97.18%, and a specificity of 78.57%. Across folds, CX3CR1 was selected in all five folds and P2RY12 in four of five folds. External validation on GSE213488 demonstrated an AUC of 0.865 (95% CI: 0.681–1.000, DeLong method; Figure 7E). Per-gene AUC values are reported in Table 3.

3.8 Drug repurposing and transcription factor prediction

Enrichr analysis identified vorinostat (SAHA), an FDA-approved HDAC inhibitor with blood-brain barrier permeability, as a hypothesis-generating drug candidate (adjusted P = 0.042; Figure 8A). This finding is consistent with preclinical studies demonstrating that HDAC inhibition restores synaptic plasticity and seizure threshold in TSC models (Basu et al., 2019). Additional candidates included piperlongumine, emetine, tamibarotene, and ciclopirox. IRF8 was predicted as the top upstream transcription factor (adjusted P = 1.45E-7; Figure 8B), consistent with its dual role as hub gene and master regulator of monocyte/microglial differentiation (Kurotaki et al., 2013).

3.9 In vitro validation in LPS-stimulated and mTOR-activated BV-2 microglia

To experimentally validate the transcriptomic findings, we examined the expression of four key hub genes in BV-2 microglial cells under two distinct activation paradigms: classical TLR4-mediated inflammation (LPS) and mTOR-driven activation (MHY1485). qPCR analysis revealed that LPS stimulation significantly upregulated Irf8 and Tlr7 mRNA, consistent with their known roles in innate immune activation (Figure 9). Notably, P2ry12 mRNA expression was significantly downregulated in the LPS group compared to Control, whereas it was maintained or slightly upregulated in the MHY1485 group. Cx3cr1 showed a decrease under LPS stimulation but an increase under MHY1485 treatment, mirroring the divergent regulation observed for P2ry12.

Western blot analysis confirmed the qPCR findings at the protein level (Figure 10). Irf8 protein was markedly upregulated in both LPS and MHY1485 groups. P2ry12 protein showed divergent regulation: decreased in the LPS group but maintained in MHY1485-treated cells. Phospho-S6 (p-S6) levels were strongly elevated in the MHY1485 group, confirming successful mTORC1 activation, but remained unchanged in the LPS group. These results provide proof-of-concept evidence that mTOR-driven activation preserves P2ry12 expression, whereas classical TLR4-mediated inflammation triggers its loss, providing mechanistic evidence that the homeostatic-to-reactive microglial transition differs fundamentally between these two activation modes.

4 Discussion

In this study, we identified a conserved transcriptomic program of innate immune activation, predominantly microglia/myeloid-associated, across FCD2a, FCD2b, and TSC. The 15 hub genes, all microglia/myeloid-associated, represent a cross-subtype disease module. Several genes warrant focused discussion.

CX3CR1 and P2RY12 are canonical homeostatic microglial markers (Hickman et al., 2013). Although often described as homeostatic markers, both CX3CR1 and P2RY12 can remain highly expressed in certain activated human microglial states, and their marked upregulation in patient tissue likely reflects an expanded microglial population rather than simple per-cell downregulation. This interpretation is supported by our in vitro experiments, which demonstrated that P2ry12 expression was preserved or upregulated under mTOR-driven activation (MHY1485) but was markedly reduced under classical TLR4-mediated inflammation (LPS; Figures 9, 10). These findings are consistent with the established role of P2RY12 in regulating microglial activation states through purinergic signaling (Haynes et al., 2006), and suggest that mTOR pathway hyperactivation, a hallmark of FCD and TSC, drives a distinct microglia/myeloid-associated program that retains features of homeostatic identity. The IRF8-IRF5-TLR7 axis represents a compelling regulatory cascade: IRF8 as master TF (Kierdorf et al., 2013), TLR7 as endogenous RNA sensor (Kawai and Akira, 2010), and IRF5 as pro-inflammatory effector (Takaoka et al., 2005). C3AR1 and complement pathway enrichment implicate complement-mediated synaptic pruning in epileptic network reorganization (Stevens et al., 2007).

Limitations include: retrospective design with public datasets; limited external validation sample size; lack of batch effect correction between training and validation datasets; and the inherent inability of bulk RNA-seq to resolve microglial abundance, activation state, or infiltrating myeloid cells. Cross-region differences between hippocampal TLE-HS and cortical FCD/TSC samples may also contribute to the observed signature, although a three-cortical-group sensitivity analysis confirmed robustness. While in vitro validation using the BV-2 murine microglial cell line provided protein-level confirmation of key hub genes, these findings require further validation in human iPSC-derived microglia or primary human surgical specimens. Additionally, drug predictions remain purely computational and require experimental confirmation. Future directions include scRNA-seq validation, IHC confirmation in surgical specimens, and functional studies using conditional knockout models.

5 Conclusion

Focal cortical dysplasia and TSC share a conserved microglia/myeloid-associated innate immune transcriptomic program. Fifteen hub genes represent a cross-subtype disease module, and CX3CR1 + P2RY12 provide a two-gene diagnostic classifier with encouraging but preliminary performance (training AUC 0.979; nested CV AUC 0.862; external validation AUC 0.865). GSEA reveals convergent TLR, NF-kB, chemokine, and complement pathway enrichment. IRF8 is the predicted master transcriptional regulator, and vorinostat is a candidate therapeutic. These findings support microglia-targeted strategies as a trans-subtype therapeutic approach in drug-resistant epilepsy.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2026.1954412/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Per-group multi-gene ROC curves. Supplementary material at the source

GO enrichment dotplots for 49 shared DEGs. Supplementary material at the source

WGCNA module eigengene dendrogram. Supplementary material at the source

Hub gene expression in GSE213488 validation cohort. Supplementary material at the source

GSEA GO BP dotplots for all comparisons. Supplementary material at the source

Cortical-only DEG overlap among FCD2a, FCD2b, and TSC (Venn diagram). Supplementary material at the source

Principal component analysis of cortical samples by disease group. Supplementary material at the source

Sample clinical characteristics (GSE256068). Supplementary material at the source

Full list of 49 shared DEGs. Supplementary material at the source

WGCNA module gene assignments (MEbrown and MEred). Supplementary material at the source

Gene significance (GS) of the 15 hub genes for TSC, FCD2a, and FCD2b. Supplementary material at the source

Complete GSEA KEGG results for the four comparisons. Supplementary material at the source

Complete GSEA GO BP results for the four comparisons. Supplementary material at the source

xCell immune scores. Supplementary material at the source

Footnotes

Edited by: Jihoon Nah, Chungbuk National University, Republic of Korea

Reviewed by: Paolo Milazzo, University of Pisa, Italy

Haifeng Ran, The Affiliated Hospital of Zunyi Medical University, China

Footnotes

https://www.ncbi.nlm.nih.gov/geo/