differential expression analysis with DESeq2 & edgeR

Question

[Rob]

Hello I have done gene differential expression analysis for HT-Seq count data using two packages DESeq2 and edgeR. I have different results forr two packages, and edgeR gives more significant genes compared to DESeq2. Why the results are different? which one is reliable? Thanks

[Kevin Blighe]

There was a question related to this just yesterday. Please confirm that you have searched for an answer? The lack of details in your post makes it very difficult for anybody to respond with any definitive answer. You need to explain exactly how you processed your data and show the code.

[Rob]

Thanks Kevin This is my code for edgeR:

#

library(edgeR)
samples<-read.csv("sampledata.csv")
expression<-read.csv("HTSeqCount.csv",header = T,row.names = 1)
###
result <- DGEList(counts=expression, genes=rownames(expression), samples=samples)
result <- calcNormFactors(result)
#
results <- result[, result$samples$Status %in% c("0", "1")]
design <- model.matrix(~as.character(Status), data=results$samples)
colnames(design) <- c("0", "1")

results <- estimateDisp(results, design)
fit <- glmFit(results, design)
lrt <- glmLRT(fit)
Final<-topTags(lrt,n=19000)
View(topTags(lrt))
#Saving results
write.csv(Final,"final.csv")

This is my code for DESeq2:

## Reading in raw data 
rdata <- read.table("HTSeqCount.txt", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
sample_org <- data.frame(row.names = colnames(rdata), c(rep("0", 22), rep("1", 22)))
colnames(sample_org) <- c("Group")
####################
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = rdata,
                              colData = sample_org,
                              design = ~Group)

dd <- DESeq(dds)
res <- results(dd)
#Saving results
write.csv(res,"res.csv")

[ATpoint]

The never ending story of edgeR vs DESeq2, this has been asked so many times before. A proper comparison requires proper knowledge of the internals, see e.g. https://mikelove.wordpress.com/2016/09/28/deseq2-or-edger/. You should also follow the manuals, as you are missing both the FilterByExpr filter for edgeR and use the relatively old LRT test rather than the currently recommended QLFTest. My honest opinion, just use what you feel more familiar with, what feels better in your hands. Both DESeq2, edgeR as well as limma-voom work well, are tested, maintained, well-accepted. I personally use edgeR because I find it a little easier to fit more complex designs via model.matrix and makeContrasts even though the same is possible with DESeq2. The latter now offers explicit fold change shrinkage which can be nice in certain situations when you really want to use the effect sizes for some ranking purposes like GSEA. Eventually, much more important is what you make out of these results the tools produce which are only the beginning of the journey.

[Rob]

Thanks The reason I am asking the difference between two packages is that when I plot heatmap with significant genes resulted from two packages I don't see pattern and I believe they are not really differentially expressed and I was thinking maybe these two packages are not appropriate tools for HT-Seq data...

[ATpoint]

So it is not about DESeq2 or edgeR but on a heatmap problem, have stated that right away. No offense but the probability is higher that there is a flaw in your heatmap code than that the arguably two most established and tested packages for DE testing produce nonsensish results.