Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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   |

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.

share|improve this question
Note the problem I mentioned in my earlier answer to your other question -- about testing only against one variable when there are other important variables (I pointed to the wikipedia article on Simpson's Paradox) - Fisher's exact test doesn't get around that. – Glen_b Dec 20 '12 at 23:59
Bioinformatics is cool!! Welcome to the site! – Kyle. Dec 21 '12 at 1:10
I'll come back and give more extensive answers soon, but the appropriate functions in R to look at will be loglin, loglm (in the package MASS, which comes with R but isn't installed by default) and glm itself. The understanding of these models will have a lot of similarity to understanding multiple regression and ANOVA - with the exceptions that the distributions aren't normal, and the logs-of-means is what the models are linear in. – Glen_b Dec 21 '12 at 22:20

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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.