Hi recently there has been some publication on the importance of GC and length bias ( Mandelboum et al, 2019, PLOS) . I'm looking into how to do this and came across packages like EDAseq. So it looks pretty straightforward with something like this.

dataOffset <- withinLaneNormalization(data,"gc",
                                      which="full",offset=TRUE)

this provide two slots, one for the normalized counts and the other for the offsets. I'm wondering if I can then use the dataOffset normalize count, say normalized_count as input for TMM normalization follow with voom. Would this work, something like.

y.df <- calcNormFactors( y.normalized_count , method = "TMM" )
voom <- voom(y.ydf )

and then just do the limma DGE?

thanks!

I did a careful examination of the phenomenon reported by Mandelboum et al (PLOS Biology 2019) and I rexamined some of our historical dataset analyses using limma and edgeR in this light. I agree that the length bias phenomenon does exist for some datasets, but it is very unlikely to cause any serious problems if you use a standard limma-voom or edgeR pipeline as described in the current documentation.

The phenomenon manifests as a expression bias for very short genes in some samples. The bias can be in either direction (up or down) for an individual sample but is consistent for all the affected genes in that sample. The latter fact induces an inter-gene correlation. The reason that limma-voom analyses are not much impacted is that:

1. the effect is too small to be important at the individual gene level and

2. at the pathway level, the limma functions roast(), fry() and camera() prevent the bias from having any impact because they correct for inter-gene correlation (as Mandelboum et al note in their paper).

Note on the other hand that preranked GSEA is sensitive to the problem and does not provide the protection that roast(), fry() and camera() do.

Regarding your actual question, I am not familiar with the output from withinLaneNormalization so I can't comment, except to say that calcNormFactors and voom expect to get actual counts and not computed quantities. As far as I can see, the EDASeq documentation does not explain what is meant by a "normalized count", but the EDASeq vignette does say that "normalized counts" can be input to edgeR and DESeq. If "normalized counts" are suitable for edgeR, then they will certainly be suitable for limma also (but I have not tried this).

Your question is a bit curious in that Mandelboum et al only finds with length bias to be a problem whereas your code deals only with GC bias and ignores length, so there is a disconnect between the results of the cited paper and your code.