Hi,

I have done differential expression on same date with both of DESeq2 and EdgeR but I am seeing contradict results. For instance DCK6 shows positive 1 fold change by DESeq2 while -1 fold change by EdgeR. I really got confused which one is right?

condition = as.factor(c(rep ("TRG1-2",26), rep("TRG4-5", 30)))
mycols = data.frame(row.names = c("A2", "A3", "A4", "A6", "A7", "A8", "A9", "A10", "A11", "A12", "B4", "B5", "B6", "B7", "B8", "B9", "B10", "B11", "B12", "C1", "C2", "C3", "C4", "C5", "C6", "G12", "D1", "D2", "D3", "D4", "D5", "D6", "D7", "D8", "D9", "D10", "D11", "D12", "E1", "E2", "E3", "E4", "E5", "E6", "E7", "E8", "E9", "E10", "E11", "E12", "F1", "F2", "F5", "F6", "F7", "G8"))
mycols$condition= c("TRG1-2","TRG1-2", "TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG1-2","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5","TRG4-5", "TRG4-5")

dds=DESeqDataSetFromMatrix(countData = df, colData = mycols, design = ~ condition) 

dds <- DESeq(dds)
DESeq2 <-results (dds, contrast=c("condition", "TRG1-2", "TRG4-5"))

> head(DESeq2)
          baseMean log2FoldChange     lfcSE     stat       pvalue      padj
CDK6     1617.5983       1.667070 0.2901847 5.744858 0.0000000092 0.0000202
CALML5    378.1458       1.496815 0.4058617 3.687992 0.0002260310 0.0653134
KLK5     1540.0631       1.472039 0.4948118 2.974946 0.0029304000 0.1895279
KRT14   22457.5507       1.337578 0.6644294 2.013122 0.0441017650 0.4903814
KRT17    9115.3139       1.182890 0.5461288 2.165954 0.0303146840 0.4473959
S100A7A   351.9248       1.049479 0.3290714 3.189213 0.0014266070 0.1329832

For EdgeR

> group = condition
> group
[1] TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2
[17] TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG1-2 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5
[33] TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5
[49] TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5 TRG4-5
Levels: TRG1-2 TRG4-5
> dim(df)
[1] 2560   56
> y <- DGEList(counts = df, group = condition) 
> y <- estimateDisp(y) 
Design matrix not provided. Switch to the classic mode.
> sqrt(y$common.dispersion)
[1] 0.6280918
> EdgeR <- exactTest(y) 
> topTags(EdgeR)
Comparison of groups:  TRG4-5-TRG1-2 
           logFC   logCPM       PValue          FDR
PPBP  -4.3340878 9.503884 3.564802e-11 9.125894e-08
CDK6  -1.5518198 8.712466 1.458599e-07 1.867006e-04
IL1B   1.7324695 9.178351 2.623373e-05 1.908504e-02
CXCL8  1.6455933 8.340310 3.129262e-05 1.908504e-02
EGR1   0.8468036 8.652308 4.432857e-05 1.908504e-02
IFIT2  0.8957873 7.535228 5.199642e-05 1.908504e-02
IL6    1.3926323 6.951407 5.218565e-05 1.908504e-02
BDNF   1.4176689 6.605966 7.471018e-05 2.134076e-02
PTGS2  1.4746062 8.352272 7.547266e-05 2.134076e-02
FOS    0.9891503 9.263358 8.336234e-05 2.134076e-02

I suspect that you've done quite different analyses in DESeq2 and edgeR, i.e., you've tested for different things in the two analyses. It's impossible to tell though, because your edgeR code is based on an object data_clean that you haven't told us anything about and which isn't used in the DESeq2 analysis.

Your edgeR code cannot be correct as you show it. It is impossible for you to have got the output you show from the code in your post because head() doesn't work on output from exactTest. You need topTags(EdgeR).

Why are you using different grouping variables for the two analyses? Why not group=condition for the edgeR analysis?

Later

Now that you have corrected your code it becomes clear that you have simply tested the same comparison in opposite directions in the two packages. In edgeR you have estimated log-fold-changes for TRG4-5 vs TRG1-2. In DESeq2 you have estimated log-fold-changes the other way around, for TRG1-2 vs TRG4-5. Naturally the logFCs will have opposite signs.