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?
expression*log(length)
)? $\endgroup$