Question: Different results with DESeq2 and EdgeR
1
gravatar for jivarajivaraj
12 weeks ago by
jivarajivaraj10 wrote:

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
edger deseq2 • 184 views
ADD COMMENTlink modified 12 weeks ago • written 12 weeks ago by jivarajivaraj10
1

I have edited your post to format your code correctly. Bioconductor is now using markdown for posts, which means that you have to insert four spaces before each line of code in order for it to be recognised as code.

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth37k

You need to show the code you used for your edgeR and DESeq2 analyses. Otherwise we have no way to tell whether either of your analyses are correct.

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth37k

Thank you I will add my code to the post, I just know the sign of top genes in two method is completely flipped :(

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by jivarajivaraj10
Answer: Different results with DESeq2 and EdgeR
1
gravatar for Gordon Smyth
12 weeks ago by
Gordon Smyth37k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth37k wrote:

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.

ADD COMMENTlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth37k

Sorry,

I have edited my post and I am sure this time I have putted same things in both tools but again CDK6 shows negative fold change :( while in DESeq2 shows positive fold change. I am seeing the sign of expression of top genes in both method is 100% different. This is HTG EdgeSeq assay with 2560 oncology biomarker genes.

> head(df[,1:5])
          A2    A3    A4    A6    A7
ACTB    8502 15629 13216 21418 32938
ATP5F1   679  1039   899  1177  1074
DDX5    8006  8063 14763 10988 10852
EEF1G  10279 15727 13000 13007 12849
GAPDH  33956 40181 49591 48761 48745
NCL    10268  8243  9375  8449  9599
ADD REPLYlink modified 12 weeks ago by Gordon Smyth37k • written 12 weeks ago by jivarajivaraj10

Again I have edited your post to display the code correctly. The Bioconductor editor requires you to indent each code line by 4 spaces.

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by Gordon Smyth37k

Thanks a lot, very helpful. Sorry as I mentioned already this is not a RNA-seq rather HTG EdgeSeq only targeting 2560 oncology biomarkers. Do you think could I use DESeq2 or edgeR for differential expression or should I use t-test (the distribution of reads is not normal though)? I saw a post mentioning voom for nanostring assay. I am not sure if I am in right path

ADD REPLYlink written 12 weeks ago by jivarajivaraj10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 109 users visited in the last hour