This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Understanding RNAseq normalization

I inherited some code from a colleague and I have doubts about whether they're doing it the right way. The input is RNAseq data, processed to get gene abundance with kallisto.

They first use the tximport package to pull in the data:

tx <- tximport(k, type="kallisto", tx2gene=tx2gene,
               countsFromAbundance="lengthScaledTPM",
               ignoreTxVersion=T);

Then they feed those counts into edgeR and normalize (convert to CPM), then run differential expression:

y <- DGEList(counts=tx$counts);
y <- calcNormFactors(y);
counts.cpm <- cpm(y, normalized.lib.sizes=TRUE)
y.classic <- estimateDisp(y)
et <- exactTest(y.classic)
. . .

From a quick inspection of the data, I think the tx$counts values are getting modified from the original kallisto esimates in that first step. So if they're already scaled, does that mean that I shouldn't be re-normalizing in the edgeR step? Or maybe that I should be using tximport differently (setting countsFromAbundance = "no"?)

It seems to me like things are getting doubly-transformed here, but I'm fuzzy on the details, even after diving into the documentation quite a bit. Help figuring out if this makes sense would be appreciated!

rnaseq normalization edger tximport
• 781 views • link
updated by ATpoint ATpoint • written by bioinfoguy2 • bioinfoguy2

1 answer

The concept of tximport is that it aggregates the transcript abundance estimates from kallisto (or salmon, stringtie, RSEM...) to the gene level. Meaning, instead of counts per transcript for a gene, you get a single count per gene. The crux is that it respects the length of the transcripts. Longer transcripts give higher counts at same expression level, and if one sample expresses longer transcripts and another shorter one, then the latter would get fewer counts, simply due to transcript composition. Tximport corrects for that, and the different options (such as lengthScaledTPM) determine how exactly this works.

In a nutshell, the code is perfectly fine and returns raw counts on the gene level that are suitable for edgeR. The edgeR code is basically fine, but the exactTest is relatively old, and generally a prefilter for low counts makes sense, check the user guide.

• 0 views • link
written by ATpoint ATpoint

Thanks, that explains a key point I was missing!

• 0 views • link
written by bioinfoguy2 • bioinfoguy2

Log in to answer this question.

More posts like this