Correct for gene length or read counts in GO enrichment analysis

Question

It is a well reported fact that GO analysis of RNAseq results is affected by a number of biases, including length bias and expression level bias.

The bioconductor package goseq allows you to correct for these biases.

By default it corrects for length bias, but you can also get it to do read count bias. Using read counts to do the correction is attractive because in theory it should account for both sources of bias ($read counts\approx expression \times length$).

I'm doing an enrichment analysis were I have tried both options (length and read counts) and get very different answers. If I run a binomial regression on expression and length vs probability of being differential, I can see that both are independently important.

> model <- glm(sig ~  expression + log(length), data=retained_genes,  family=binomial(link="logit"))
> print(anova(model, test="Chisq"))

Analysis of Deviance Table

Model: binomial, link: logit

Response: sig

Terms added sequentially (first to last)

                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                                    6676     4507.1              
expression              1  114.998      6675     4392.1 < 2.2e-16 ***
log(length)             1  102.553      6674     4289.5 < 2.2e-16 ***
expression:log(length)  1   34.094      6673     4255.4 5.252e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Answer 1

My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and before enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.

Now, to the question itself:

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.

Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.

Answer 2

I am not sure if it makes sense to use read counts as bias instead of gene length (and I certainly wouldn't expect the same results).

Do you use total read counts of all your samples (library size)?

The correction for gene length is a pure technical one, the longer a gene the more reads will align (and higher read count genes are easier significant, since they are way above the noise threshold). If you use read counts, you also have a biological factor (expression) in there, which (I think) is the stuff you test with statistics (e.g., with edgeR) and thus not the bias you want to correct for.