# Correct for gene length or read counts in GO enrichment analysis

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

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?

• Did you test for the interaction of them (expression*log(length))?
– llrs
Aug 21 '17 at 15:37
• I had not. Now I have. See updated question. Aug 21 '17 at 15:51

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.

• These bias' work through power. You cannot magically generate power out of thin air, so there are two ways you could remove the bias pre or at the DE stage: 1) reducing variance between samples 2) down-weighting more highly powered genes. As Mike Love notes in that link, modeling in the row-wise read counts is not compatible with DESeq2s methods (it were possible to adjust for extra power of highly expressed genes at the DE stage, presumably DESeq would already do it). Also who does one measure the accuracy of DEG prediction with no ground truth? Aug 21 '17 at 16:04
• Taking into account the length of the genes won't reduce the variance between samples, so your only option is to down-weight highly powered genes. Maybe there are other methods to take it into account (or you could develop them). To find the true DEG for a given comparison: You can explore the FDR values of the genes in the comparison. The FDR values are assumed to follow an uniform distribution, because each gene has the same probability to be DEG. By exploring how uniform is that distribution you can find if the model captures the biases or not. You could also use the SeqGSA method
– llrs
Aug 21 '17 at 18:56
• But I'm not interested in finding a better set of DEGs, rather, in better assessing whether their division in categories is what you would expect by chance. GOSeq allows you to do this while correcting for bias, but only one bias at once. Aug 22 '17 at 8:45
• The problem in this approach is that the DEG you test, are affected by this length and GC bias. If you want to to assess their division in categories use topGO, because it is the only method which does an ORA analysis taking into account the structure of the GO graph. GOSeq is biased by not taking into account the relationships between GO terms.
– llrs
Aug 22 '17 at 8:56

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.

• Imagine you have two genes: Gene A has 10,000 counts in Condition A and 5,000 in condition B. Gene B has 1,000 counts in Condition A and 500 in condition B. The lfcs are the same, but Gene A is more likely to be call differential than gene B. Now consider two GO categories. Category X is made up of 10,000/5,000 genes and Category Y is made up of 1,000/500 genes. Cat X is more likely to be enriched than Cat Y, despite the real LFC in the genes being the same. Aug 21 '17 at 13:27
• I don't know if the second example is true (that category X is more likely enriched). Maybe you should ask real statisticians on bioconductor?
– benn
Aug 21 '17 at 14:09