Question: Unexpected outcome: Gene Enrichment using Fisher's Test
4
gravatar for Adrian Pelin
4.2 years ago by
Adrian Pelin2.2k
Canada
Adrian Pelin2.2k wrote:

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.

 

 

ADD COMMENTlink modified 4.2 years ago by aodaaodaslman20330 • written 4.2 years ago by Adrian Pelin2.2k
3
gravatar for Devon Ryan
4.2 years ago by
Devon Ryan86k
Freiburg, Germany
Devon Ryan86k wrote:
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 :)

ADD COMMENTlink modified 4.2 years ago • written 4.2 years ago by Devon Ryan86k
3
gravatar for russhh
4.2 years ago by
russhh4.1k
UK, U. Glasgow
russhh4.1k wrote:

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.

ADD COMMENTlink modified 4.2 years ago • written 4.2 years ago by russhh4.1k

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.

ADD REPLYlink written 4.2 years ago by Devon Ryan86k
2
gravatar for Michael Dondrup
4.2 years ago by
Bergen, Norway
Michael Dondrup45k wrote:

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 

 

ADD COMMENTlink written 4.2 years ago by Michael Dondrup45k

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?

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by Adrian Pelin2.2k
1

no it doesn't include FDR correction

ADD REPLYlink written 4.2 years ago by russhh4.1k

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.

ADD REPLYlink written 4.2 years ago by Adrian Pelin2.2k
2

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 micorarray probe mapping one-mahy to multiple ribosomal genes.

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by russhh4.1k

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?

ADD REPLYlink written 4.2 years ago by Adrian Pelin2.2k
1
p.adjust(p, method = p.adjust.methods, n = length(p))
ADD REPLYlink written 4.2 years ago by russhh4.1k

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.

ADD REPLYlink written 4.2 years ago by Adrian Pelin2.2k
2

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

ADD REPLYlink written 4.2 years ago by russhh4.1k

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

ADD REPLYlink written 4.2 years ago by Adrian Pelin2.2k

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.

ADD REPLYlink written 4.2 years ago by karl.stamm3.4k
1

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.

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by Devon Ryan86k
1
gravatar for EagleEye
4.2 years ago by
EagleEye5.9k
Sweden
EagleEye5.9k wrote:

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.

ADD COMMENTlink written 4.2 years ago by EagleEye5.9k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1280 users visited in the last hour