Unexpected outcome: Gene Enrichment using Fisher's Test

Question

[Adrian Pelin]

Hello,

I have an enrichment analysis I am trying to work out and it gives an unexpected result and I was wondering if I am doing it correct.

I have found KEGG categories for all genes (Ref. Set) in the genome of my organism, when I add up the number of genes in each category, the total is 3426.

I then have chosen a subset of these genes (Test Set) based on a criterion of my choice, and would like to see if these are enriched for a specific pathway/function. Finding KEGG categories for these results in a total of 98 categories.

I then notice that the biggest category in my "Test Set" is "Ribosome", with 12 genes (out of 98) being part of that category. When I look at my "Ref. Set", I see that overall in the genome 87 (out of 3426) are part of that category.

Just looking at ratios:

12/98 = 0.1224
87/3426 = 0.02539

It looks like there is more of that category in the test set compared to reference set. Now to the Fisher's Test:

> fisher.test(matrix(c(12, 75, 86, 3253), 2, 2), alternative='less')$p.value
[1] 0.9999993

It looks like the p-value is very high, indicating no enrichment whatsoever. What am I missing here? Thank you.

Answer 1

[Devon Ryan]

fisher.test(matrix(c(12, 75, 86, 3253), 2, 2, alternative="greater")$p.value

Given how your matrix is set up you just need to look at the other tail :)

Answer 2

[russhh]

You're looking for enrichment within your gene set, you should use alternative = 'greater'. This tests whether the odds ratio for annotation of your gene set to the ribosome is greater than it would be under the null hypothesis (no association): effectively, are there more ribosomal genes present than expected.

A simple sanity check that you can do is to compare the results of the Fisher test when an extra one of your genes is annotated to the ribosome:

fisher.test(matrix(c(12, 75, 86, 3253), 2, 2), alternative='less')
# pvalue = 1; OR = 6.044871
fisher.test(matrix(c(13, 74, 85, 3254), 2, 2), alternative='less')
# p-value=1; OR = 6.715563

An additional ribosomal gene increased the odds ratio but had no effect on the pvalue.

Try it with alternative=greater:

fisher.test(matrix(c(12, 75, 86, 3253), 2, 2), alternative='greater')
# p-value = 4.573e-06; OR = 6.04
fisher.test(matrix(c(13, 74, 85, 3254), 2, 2), alternative='greater')
# p-value = 6.86e-07; OR = 6.715

Using the correct test, an additional ribosomal gene increased the odds ratio and also decreased the pvalue.

[Devon Ryan]

I second the idea of a simple sanity check. It's easy enough to not match the matrix orientation that fisher.test() is expecting that you're pretty likely to muck things up otherwise.

Answer 3

[Michael]

You have to use alternative="greater", not "less", because the odds ratio is calculated as odds(1.column)/odds(2.column). And for your case, your test-set is in the first column and you want to test whether the odds to see a category X in test is significantly greater than in all genes, thus whether the odds ratio is > 1. Your odds are approximately:

(12/75) / (87/3253)
[1] 5.982529

That's similar to the estimated odds-ratio from the test:

> fisher.test(matrix(c(12, 75, 86, 3253), 2, 2), alternative='greater')

    Fisher's Exact Test for Count Data

data:  matrix(c(12, 75, 86, 3253), 2, 2)
p-value = 4.573e-06
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
 3.245896      Inf
sample estimates:
odds ratio 
  6.044871

[Adrian Pelin]

Oh, sounds simple, I suppose the last parameter is very important:)

Thanks for the help Devon Ryan, Michael and russ_hyde, appreciate it!

Edit: I was wondering, does this test include FDR correction?

[russhh]

no it doesn't include FDR correction

[Adrian Pelin]

What purpose would FDR testing serve in my case? When I do the same analysis in blast2go based on GO terms, I find the same category enriched with a p-value 3.1e-4, which is great, but if I include FDR I see no category enriched. I am wondering what could be the reason and how is FDR done in this case. If I try to present this data to biologists, would I end up telling them that it's enriched for Ribosomal genes or would I say it is not.

[russhh]

I imagine there are many more GO terms considered in that program than the 98 KEGG pathways you considered here, so a lot of genuine associations may be stamped out by the FDR correction within blast2go. Anyway, you should treat FDR results with kid gloves when there are high levels of dependence between different pathways, something that is inevitable with GO, given it's hierarchical structure. Perhaps you could present before and after correction values for an interpretable number of your best-scoring pathways.

Here, your ribosome result will be resistant to FDR testing (it's modified pvalue will be between 10^-3 and 10^-4) if you adjust for the 98 pathways. However, you've just given away that you tested more pathways than that, so perhaps you should FDR correct for 98 + the number of GO terms you also tested....

I think it's more important that you check that the pathways that appear to be fisher-significant have come from technically independent assays: I mean, it'd be pretty bad if you had many ribosomal genes in your test set due to a single microarray probe mapping one-many to multiple ribosomal genes.

[Adrian Pelin]

Thanks, that's informative. I am looking at a subset of genes with no heterozygosity present in them, data collected from a pop-gen analysis. How would I go by manually correcting p-values with FDR testing? Is there an R command?

[russhh]

p.adjust(p, method = p.adjust.methods, n = length(p))

[Adrian Pelin]

Alright so this is actually pretty simple to do, I think I understand it better. What I do is I test every single category that has at least 1 gene enriched for it from my "Test Set".

I think it does not make sense to include all p-values in this correction, since I would discard some even if they were significant. For instance, one gene in one category being present in my Test Set could show a significant p-value if there is a total of only 2 genes in those categories in the Reference Set. I think I should set the cutoff at at least 2 genes in a pathway before I calculate all p-values, then do correction on that.

From what I understand BH and FDR are identical, and are recommended, so I will stick to that. One more thing I was wondering about. I can enrich with KEGG, GO and KOG. KOG seems to have the least number of categories, and it bins genes with more than one pathway into "Multiple Classes", so the amount of genes in one category such as "Transcription" does not overlap genes in other Categories. It seems that based on this, it would be the most reliable way of testing for enrichment since there is minimal nr of categories and no overlap.

[russhh]

You have to be careful when you analyse a dataset after having seen the results... ;)

[Adrian Pelin]

Agreed. I appreciate your help and the help of others a lot, thanks for guiding me through this!

[karl.stamm]

Minimal # of categories and no overlap sounds like a benefit for the algorithm but a detriment to your final objective.

Redundant annotations can make more interesting categories. We know TFs are used in dozens of places, as are structural proteins. So when it says "Kidney Nephron" as a GO Term, the process generalizes to other unknown areas. I couldn't do anything with a result that said "Transcription" is impacted by Drug X. It's too vague.

[Devon Ryan]

If you're going to test all of (or even a decent portion of) the categories then you'll need to adjust the resulting p-values. For how that works, it depends on the algorithm you use. Generally one would use BH correction, which is probably recommended.

Edit: For why we need to adjust p-values, this XKCD comic is surprisingly good.

[seta]

Hi Michael,

I used the similar fisher test on SNPs (2400 SNPs), but when I calculated OR with fisher test, there are multiple SNPs with OR < 1 while "alternative=greater" was used in the fisher test command; so I expected all enriched SNP (OR > 1). Could you please let me know what's happened?

Thanks

Answer 4

[EagleEye]

You can have a look at this post if you have genes from Human :

Gene Set Clustering based on Functional annotation (GeneSCF)

This will help you for sure.