DESeq2 with different replicates

Question

[robertorun]

There are 104 samples, 38 samples are treated, 66 samples are control. The following script was used to do DGE analysis. But there is no any different expression genes, even though using the argument of dds$padj<0.05, alpha = 0.05, lfcThreshold=1. But there are more than 6000 DGEs, if you use edgeR, DEGSeq or GFOLD et al.

Would you like to tell me where is wrong? Or which soft package is more fit to do this analysis (with different replicates)? Thanks very much!

[robertorun]

code is the following:

library('DESeq2')
directory <-"./At_Count/"
sampleFiles <- grep("*.txt",list.files(directory),value=TRUE)
sampleFiles

sampleCondition <- c("A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.I","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W","A.W")
sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition)
sampleTable

ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~condition)
ddsHTSeq

colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c("A.I","A.W"))
#ddsHTSeq <- ddsHTSeq[rowSums(counts(ddsHTSeq)) >520,]
ddsHTSeq

dds<-DESeq(ddsHTSeq)
summary(dds)
table(dds$padj<0.001)

res<-results(dds,alpha = 0.05, lfcThreshold=2)

summary(res)
table(res$padj<0.05)

write.csv(as.data.frame(res),file="./At.IW.deseq2.MAplot.LFC2.Alp05.csv")

png(file="./At.IW.deseq2.MAplot.LFC2.Alp05.png",width=10,height=7.5,units="in",res=600)

layout(matrix(c(1,2,3,4,5,6), 1, 1, byrow=TRUE))
plotMA(res,ylim=c(-10,10),main="Improved VS Wild of At")

dev.off()

[pld]

Maybe drop the lfcThreshold argument from results. When this is supplied, DESeq2 performs a test to determine if the absolute value of the fold change is greater than the value of lfcThreshold. This is totally different than asking if a gene is simply differentially expressed.

When lfcThreshold is supplied, the null hypothesis is "between these two conditions the absolute value of the fold change of a gene is not greater than lfcThreshold". When you're just looking looking at p-values, your null hypothesis is "gene x is not differentially expressed."

Just drop it out and filter the results on log2FC after.