I have some RNASeq data, that has been spiked-in with exogenous ERCC controls. I have the counts from htseq-count and now I want to normalize + test for differential expression. For this, I am using two methods: DESeq and limma+voom. Below is my code for both the tests. Can someone tell me if what I am doing is correct or not? I do think that I am slightly mistaken in the voom code.
##### Using DESeq library(DESeq) # get the size factors from the spiked-ins # dat.spike.in is count data of only spike ins cds <- newCountDataSet(countData = dat.spike.in, conditions = as.character(condition)) cds <- estimateSizeFactors(cds) nf <- sizeFactors(cds) # add the size factors to data without spike-ins # dat.without.spike.in is count data where spike-ins are removed cds <- newCountDataSet(countData = dat.without.spike.in, conditions = as.character(condition)) sizeFactors(cds) <- nf # differential expression test cds <- estimateDispersions( cds, method= 'pooled', sharingMode='maximum', fitType='local') deseq.res <- nbinomTest( cds, "A", "B") ##### Using Limma + Voom library(limma) library(edgeR) groups = as.factor(condition) # conditions design <- model.matrix(~0+groups) colnames(design) = levels(groups) nf <- calcNormFactors(dat.spike.in) # calculation norm factors using spike ins only voom.norm <- voom(dat.without.spike.in, design, lib.size = colSums(dat.spike.in) * nf) # add that to voom
Awaiting your input and suggestions!
P.S: Crossposted on Biostars to get inputs from both communities.