Gene set enrichment analysis on differential phosphorylation sites

Question

I have:

• A list of differentially phosphorylated sites in a knockout condition. Some genes contain as many as 70 possible phosphorylation sites; others contain only one.
• A list of genes belonging to a specific gene set annotation.

How can I test the differentially phosphorylated proteins for enrichment of this annotation?

A few ideas I’ve considered:

1. Ignore the number of phosphorylation events detected within a gene and simply count the gene as differentially phosphorylated if it contains at least a single site that is differentially phosphorylated. Compare the enrichment score for this selected set to the enrichment score for the set genes containing at least one site that is not differentially phosphorylated. The problem here is that genes with a large number of phosphorylation sites have almost no influence in the enrichment score, since it’s almost certain that they’ll have at least one site that is differentially phosphorylated and at least one that is not.
2. Mark each candidate phosphorylation site as either “in the set” or “not in the set” based on the protein in which it’s found. Then perform the enrichment analysis using the set of annotated phosphorylation sites instead of using the traditional enrichment analysis performed at the gene level. The potential problem with this approach is that it may place too much influence on genes with many potential phosphorylation sites.
3. Aggregate all candidate phosphorylation sites within a gene and use some numerical threshold to determine whether the gene is differentially phosphorylated or not. (There are various ways that this could be done.) Then perform enrichment analysis using the resulting set of differentially phosphorylated genes. A possible problem here is that some of the phosphorylation sites may be more functionally important than others, so it’s not clear how to weight the relative importance of individual phosphorylation sites within a gene.

I realize the goal here isn’t well-defined mathematically; I’m mainly curious what approach makes the most sense given the biological context. Of the above approaches, I’m currently leaning towards approach 2 because it’s straightforward to implement and it at least attempts to account for the variable number of phosphorylation events within a gene.

NOTE: I also have normalized phosphoprotein and protein abundances for all of these sites, obtained from mass spec. So if the solution requires an alternative method of computing differential phosphorylation, that’s fine.

Answer 1

Your dataset looks perfect for the SUMSTAT enrichment test.

1. 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.
2. 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.
3. 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:

1. 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.
2. 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).