Enrichment analysis by gene duplication level

Question

Biological Background

Over time, some plant species tend to duplicate their entire genomes, gaining an additional copy of each gene. Due to the instability of this setup, many of these genes are then deleted, and the genome rearranges itself and stabilises, ready to duplicate again. These duplication events are associated with speciation and invasion events, and the theory is that the duplication helps plants to adapt faster to their new environments.

Lupinus, a genus of flowering plant, invaded the Andes in one of the most rapid speciation events ever detected, and what's more, it seems to have more duplicate copies in its genome than the most closely related genus, Baptisia.

And now the mathematical problem:

The genomes of a member of Lupinus and a member of Baptisia have been sequenced, providing raw data about 25,000 genes in each species. By querying against a database of genes of known function, I now have a "best guess" for what functions that gene might perform - so for example, Gene1298 might be associated with "fructose metabolism, salt stress response, cold stress response". I want to know, if there was a duplication event between Baptisia and Lupinus, whether gene loss took place at random, or whether genes performing particular functions were more likely to be kept or deleted.

I have a script which will output a table like the one shown below. L * is a count of all Lupinus genes associated with the function. L 1+ is a count of lupinus genes associated with the function where at least one duplicate copy exists. I can get it to produce L 2+, L 3+ etc., although L 1+ is a much more reliable group than L 2+ due to the sequencing process.

Function            | L *  | L 1+ | B *  | B 1+ |
fructose metabolism | 1000 | 994  | 1290 | 876  |
salt stress         | 56   | 45   | 90   | 54   |
etc.

What I would like to do is to test, for each gene function, whether there are more or fewer genes with duplicates than might be expected purely by chance in Lupinus and Baptisia, and whether Lupinus differs from Baptisia in the ratio of observed to expected.

The best thing I have so far

Previous studies on different species have used Enrichment Analysis, with Fisher's Exact Test and FDR correction for multiple sampling, to do a contingency test on each row.

It would be nice to improve on this; I'm not sure this sounds like the best way to do it.

Glen_b has suggested using a GLM to analyse the data; I have played around with GLMs in JMP8, which has been interesting, but I will admit to not really understanding them.

That said, I'm trying to use R instead now.

What am I using this for?

This was originally supposed to be as part of a short research project I'm doing at university, but has now spanned off into an enormous genome annotation project. Why? Because bioinformatics is cool. Being able to take a string of A,T,C and G and use it to infer information about events which happened millions of years ago is amazing.

Needless to say, I am not going to try and submit any kindly provided answer as my own work. I would be happy to include an acknowledgement in the paper if I use a method suggested here in the submitted work.

Answer 1

While I agree that Fisher's test (or something similar) may be the most natural approach here, how about this:

1. For each unique gene, you determine the difference in number of duplications in L and B
2. Order genes by this difference. Now the genes showing most differences between the species will be on the top of your list.
3. Apply a gene set enrichment test to the ordered list of genes. For example, you can use a modified Fisher's method from my package tmod, for which you would have to define your gene sets (it should be quite straightforward). Note that Fisher's method is not related to Fisher's test.

The modified Fisher's test (dubbed CERNO by the authors who first described it in this context) takes any ordered list of genes as the input, as long as you can group them in some useful categories.

The advantage of this approach is that apart from a p-value, you can easily calculate the effect size of the enrichment and visualize it (for example, as a ROC curve over the ordered list of genes). This gives you a much better idea how much what you observe really matters for the biology you study.

Answer 2

As you say, you're asking two distinct questions.

Question 1 "is the ratio L*/L1+ different from B*/B1+ for a given gene function"

this might be best answered with Fisher's exact test using data across the row as you found previously.

Question 2 "is the ratio: genes where there is a single copy / genes where there is more than one copy, different between gene functions?"

I think this also might be best answered with Fisher's exact test. You'd test the ratio of L*/L1+ for gene function 1 against L*/L1+ for gene function 2. Then gene function 1 Vs gene function 3, etc.

Neither of these sets of questions gets at whether or not they are being maintained/deleted faster than expected purely by chance, only whether they are being deleted/maintained at rates different from each other. To know if they were being deleted/maintained at a rate different than by chance, you'd need to know the ratio of singlecopy/multiplecopies for lots of DNA regions that are only being affected by chance. If you could find such regions, you'd end up with a "Function group" where function is "None". You'd then compare this to your other gene function groups in the same way as I described in question 2.