Gene set enrichment analysis on differential phosphorylation sites
Your dataset looks perfect for the SUMSTAT enrichment test.


You need to come up with a statistic representing your gene. The simplest ideas here are the number of sites or, probably better, a proportion of phosphorylated sites.
Now you can compute a statistic for every gene set, for example just sum (SUMSTAT). You can also have average, sum of squares or something else.
Geta a null distribution of your gene set statistics by permutations. You can either randomize statistics for genes, or just randomly assign phosphorylated sites across the genome keeping the total number of sites.


Now you can compute p-value just by comparing your value with the null-distribution.

You have to be cautious of two things:


You are testing multiple gene sets, i.e. you perform a statistical test for every gene set. Therefore correct your significance for multiple testing. I suggest using FDR for this.
Beware of potential biases when doing permutations. The obvious one is gene length, i.e. longer genes have higher chance of getting at least one site. But there can be other ones like GC-content, chromosomes, etc. You can overcome this by using more realistic permutations, or by controlling for the potential correlations. You can get a few ideas on fighting biases from this paper (sorry for the self-advertisement).