I have gene expression data from RNA from isolated brain tissue with multiple experimental groups (control and treatment, males and females, and two brain regions).

I already have lists of genes ranked based on differential expression between the control and treated groups with p-values and log fold changes (via DEseq).

I would like to conduct gene set enrichment analysis (GSEA), then create a single graph that compares the magnitude and significance of the most highly enriched terms for each the four groups (sex and brain region) between each other. I am imagining a single dotplot with GO terms on the y-axis and the four groups on the x axis.

I can't find an R package or software that can do this- is an analysis like this possible and if so what are the steps?