How do I calculate the p-value in a gene enrichment analysis, using parametric and non-parametric methods?

Question

I have a list of 20 QTL from an experiment involving a certain phenotype related to bone development. Of all the candidate genes within the 20 QTL intervals, 5 are associated with a certain disease in humans. I tested the set of 20 QTL to see if they are enriched for the disease-associating genes. Each of the 20 QTL has a certain length in base pairs. I took that set of 20 lengths, and randomly placed them in the mouse genome. I forced them not to overlap and to land no closer than 500,000 bases on either side. Then I took note of the genes inside each interval and counted how many matched the names contained in the disease-associating list of genes. I repeated the process 10,000 times. This provides the expectations given randomness... how many times can I expect to see disease-associating genes inside a set of 20 randomly placed QTL? I compare that to my actual observed number. Then I calculate the probability.

To do that, I simply counted how many times 1 disease-associating gene appeared within each set of 20 randomly placed QTL and divide by 10,000. That gives me the frequency. I did that for 2 disease-associating genes, and 3, and so on up to 5, which is my observed value. I find that the frequency of 5's in the random data is less than 0.05. Is this appropriate?

Also, I see online some might suggest a parametric method, like Fishers' exact test, but that is supposed to be used for very small data sets, am I correct? Am I missing something there? How can I implement a parametric calculation for p?

Answer 1

The procedure you described has been implemented before in a package called GAT (Genomic Association Test), it also takes account of various things that might confound such an analysis, the presence of isochores in the genome, or gaps in the genome assembly.

There are parametric alternatives, but they don't account for such bias'. Fisher's exact test is only used on small data sets because it is computationally intensive, rather than because its not valid. The hypergeometric test provides a less intensive way to do the test and should produce the same answer.

The hypergeometric test is based on drawing balls from a urn. The balls are two colours and the test asks if I draw n balls, what is the probability that x of that will be red. Here balls represent bases in the genome, red ones are bases in disease associated genes and n is the number of bases in your QTLs. The only downside of this is that addition to not controlling for isochores, genomic gaps etc, it will be confounded by clustering of functionally related genes in the genome, and possibly by gene length effects. Try R's phyper function.