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.