Sorry, I was unclear in my question.
I know the DESeq2 vignette but in the EDAseq I didn't find how to handle DESeq2 while there is EgdeR and DESeq.
I was wondering if there was something like this (reported in EDAseq):
## ----deseq------------------------------------------------
library(DESeq)
counts <- as(dataNorm, "CountDataSet")
sizeFactors(counts) <- rep(1,4)
counts <- estimateDispersions(counts)
res <- nbinomTest(counts, "wt", "mut")
head(res)
BUT for DESseq2.
What I do is:
## ----deseq2------------------------------------------------
library(DESeq2)
# data is data from EDAseq example
dataOffset <- withinLaneNormalization(data,"gc",
which="full",offset=TRUE)
dataOffset <- betweenLaneNormalization(dataOffset,
which="full",offset=TRUE)
dds <- DESeqDataSetFromMatrix(countData = dataOffset@assayData$normalizedCounts,
colData = pData(dataOffset),
design = ~ conditions)
EDASeqNormFactors <- exp(-1 * offst(dataOffset))
normalizationFactors(dds) <- EDASeqNormFactors
# dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds)
res <- res[order(res$pvalue),]
head(res)
Thank you,
P
Hi Pietro, please use the "comment" function to comment on an answer, and the "answer" function only to add an answer to the original question.
I see what you mean now: you are looking for a coercion method to
DESeqDataSet
similar to the one that we have forCountDataSet
.There isn't one because nobody ever asked, but I can certainly add it to the
EDASeq
package.Note that the
DESeq
section in theEDASeq
package is a workaround due to the fact thatDESeq
could not handle gene-specific offsets. SinceDESeq2
can, the recommended pipeline is now similar to the one described in theedgeR
section.Namely, computing offsets with
EDASeq
and pass them toDESeq2
, with no need to use the normalized counts.In your example, you need to modify the code to use the raw counts to create the
DESeqDataSet
otherwise, you are normalizing twice. i.e., changing the line that creates thedds
object to:Thanks a lot!
Everything much clear now!
Pietro
Sorry to bother but I'd like to know if it's correct to pass the normCounts(dataOffset) obtained in the previous example directly to voom in order to to process RNA-Seq data prior to linear modelling in limma.
I'm trying to make a comparison on the different pipelines to make a DEA from my RNAseq data.
Thanks in advance,
P
## ----voom-limma-----------------------------------------------
logCPM <- limma::voom(normCounts(dataOffset), design)
cont.matrix <- limma::makeContrasts(contrasts = contr, levels = design)
fit <- limma::lmFit(logCPM, design)
fit <- contrasts.fit(fit, cont.matrix)
fit <- limma::eBayes(fit, trend = FALSE)
DEGs <- limma::toptable(fit, coef = 1, adjust.method = "fdr",
number = nrow(normCounts(dataOffset)))
It looks reasonable, but we haven't tested this approach, so I would carefully look at the results, residuals, etc. to make sure that everything looks reasonable in terms of model fit. It might be good to check with Gordon Smyth or other limma developers.
If you have a question about limma, I'd make a new post. But please note the advice here:
A: edgeR and DESeq2