8
$\begingroup$

I have RNA seq data which I've quantified using Kallisto. I'd like to use tximport to transform the read count data into input for EdgeR, following the R code supplied in the tximport documentation:

cts <- txi$counts
normMat <- txi$length
normMat <- normMat/exp(rowMeans(log(normMat)))
library(edgeR)
o <- log(calcNormFactors(cts/normMat)) + log(colSums(cts/normMat))
y <- DGEList(cts)
y$offset <- t(t(log(normMat)) + o)
# y is now ready for estimate dispersion functions see edgeR User's Guide

I would ideally follow this up with EdgeR's exactTest. However, I'd also like to remove genes with low read counts using code such as this:

filterCounts <- function(counts) {
  cpms <- cpm(counts)
  keep <- rowSums(cpms > 1) >= 1
  counts[keep,]
}

The only issue is that I can't tell whether I should remove the genes with low read counts before or after I normality the counts while importing with tximport in the first code. I'm thinking that it would be best to do this before, editing out the abundance, counts, and length rows of the txi data frame which correspond to lowly expressed genes, but I would like this to be confirmed before I proceed.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

The more genes you have the more robust the scaling factor is (among the reason why one doesn't normalize to ERCC spike-ins without a compelling reason), so I suppose in theory it's better to filter after determining the scale factor. Having said that, I'd be surprised if the results changed much either way. Unless you end up filtering out a LOT of genes the scaling factors should be fairly robust to the change.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Since edgeR uses a robust estimator for the scaling factor (TMM?), you’re right, it shouldn’t make much difference. $\endgroup$
    – Konrad Rudolph
    Commented Jul 27, 2017 at 14:13
  • $\begingroup$ tximport's documentation implies that I need to supply my own offset when moving data into edgeR. If I create an offset matrix using my above code and then remove a large portion of lowly expressed genes, will the previously created offset matrix be incorrect? Do I even need to supply an offset in this case as I think edgeR can calculate it automatically. $\endgroup$
    – J0HN_TIT0R
    Commented Jul 27, 2017 at 15:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.