Different results with DESeq2 and EdgeR
1
1
Entering edit mode
AZ ▴ 30
@fereshteh-15803
Last seen 11 months ago
United Kingdom

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
desEQ2 edger • 1.2k views
ADD COMMENT
1
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
1
Entering edit mode
@gordon-smyth
Last seen 19 minutes ago
WEHI, Melbourne, Australia

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 COMMENT
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 739 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6