Search
Question: DESEQ2: apparently wrong log2(fold-change) calculation
0
gravatar for Osvaldo
19 months ago by
Osvaldo20
Spain
Osvaldo20 wrote:

hi there!

My question is related to a DESeq2 output where apparently log2 fold-changes are wrongly calculated.

 

Let me, please, explain the case with one gene result:

Raw counts for that gene (HTSeq-count):

Ctrl   2200, Treatment 2097

Normalised counts by DESeq2:

Ctrl 2215.3490440231           Treatment 2082.4709372308

Log2(treatment/control) from differential expression test

Name     baseMean     log2FoldChange      lfcSE      stat     pvalue        padj
AAAS     2148.9099906269       11.0693893399             0.0568426632           194.7373455221            0              0

I expected the log2fc to be -0.0892376619

I must say that this test is done without replicates, only one control sample vs its corresponding treament sample. Could this affect in some way??

When I perform the same test with 3 controls and 3 treatments, used as replicates for the two conditions, the log2(fc) are correctly calculated.

**The control and treatment samples shown here participate also as one of these 3 replicates for each group.

 

thanks very much !

 

 

 

ADD COMMENTlink modified 19 months ago • written 19 months ago by Osvaldo20
0
gravatar for Michael Love
19 months ago by
Michael Love14k
United States
Michael Love14k wrote:
Can you show your code? Note that 2^11 = 2048, so I'd guess you didn't supply a design formula and are testing the intercept. This information would also show at the top of the results table when printed in R.
ADD COMMENTlink written 19 months ago by Michael Love14k
0
gravatar for Osvaldo
19 months ago by
Osvaldo20
Spain
Osvaldo20 wrote:

Sure I can show the code. I am pasting it below for this case. Thanks very much !

 

library("BiocParallel")
register(MulticoreParam(2))
library("DESeq2")

countTable = read.table("Treatment_3_vs_Ctrl_3.txt",header=TRUE, row.names=1)
condition = factor(c("Ctrl_3","Treatment_3"))
libType = c("oneFactor","oneFactor")
condition <- relevel(condition, "Ctrl_3")
experiment_design=data.frame(row.names = colnames(countTable), condition, libType)


cds <- DESeqDataSetFromMatrix(countData = countTable, colData=experiment_design, design=~condition)
cds_DESeqED <- DESeq(cds,parallel = TRUE)
res <- results(cds_DESeqED,parallel = TRUE,alpha = 0.05, pAdjustMethod = "BH")
write.table(res,file = "Treatment_3_vs_Ctrl_3.differentialExpression.txt",row.names = TRUE,col.names = NA,append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".")


normalizedReadCounts = counts(cds_DESeqED,normalized=TRUE)
write.table(normalizedReadCounts,file = "Treatment_3_vs_Ctrl_3.normalizedCounts.xls",row.names = TRUE,col.names = NA,append = FALSE, quote = FALSE, sep = "\t",eol = "\n", na = "NA", dec = ".")

ADD COMMENTlink written 19 months ago by Osvaldo20
What is printed above res when you print to console?
ADD REPLYlink written 19 months ago by Michael Love14k
0
gravatar for Osvaldo
19 months ago by
Osvaldo20
Spain
Osvaldo20 wrote:

please find below what is printed out. Thanks!

 

cds_DESeqED <- DESeq(cds,parallel = TRUE)
estimating size factors
estimating dispersions
gene-wise dispersion estimates: 2 workers
mean-dispersion relationship
final dispersion estimates, MLE betas: 2 workers
fitting model and testing: 2 workers
Warning message:
In checkForExperimentalReplicates(object, modelMatrix) :
  same number of samples and coefficients to fit,
  estimating dispersion by treating samples as replicates.
  read the ?DESeq section on 'Experiments without replicates'
> res <- results(cds_DESeqED,parallel = TRUE,alpha = 0.05, pAdjustMethod = "BH")

 

ADD COMMENTlink written 19 months ago by Osvaldo20

I mean, what is printed when you type:

> res [enter]

Basically, I suspect that you may have ran DESeq() with design=~1 to get the previous result.

ADD REPLYlink written 19 months ago by Michael Love14k

hi,

I've figured it out, there is a bug in the design replacement in the parallel routine in DESeq(). I'll fix this in the devel branch, but in the mean time, you can just not use the parallel argument to avoid the bug.

ADD REPLYlink written 19 months ago by Michael Love14k

Fixed in v1.11.30

ADD REPLYlink written 19 months ago by Michael Love14k
0
gravatar for Osvaldo
19 months ago by
Osvaldo20
Spain
Osvaldo20 wrote:

yep, sorry.... here you go:

 

> res
log2 fold change (MAP): Intercept
Wald test p-value: Intercept
DataFrame with 23368 rows and 6 columns
              baseMean log2FoldChange      lfcSE       stat        pvalue
             <numeric>      <numeric>  <numeric>  <numeric>     <numeric>
A1BG          5.496657       2.458554 0.97249497  2.5280895  1.146851e-02
A1BG-AS1      2.000048       1.000035 1.58919037  0.6292732  5.291702e-01
A1CF        172.955513       7.434257 0.17722342 41.9485028  0.000000e+00
A2LD1       139.520753       7.124336 0.19661333 36.2352646 1.695768e-287
A2M         647.688886       9.339157 0.09746897 95.8167192  0.000000e+00
...                ...            ...        ...        ...           ...
ZYG11B      1126.95771      10.138218 0.07277373  139.31150  0.000000e+00
ZYX         1775.77177      10.794230 0.06014082  179.48261  0.000000e+00
ZZEF1       2234.31414      11.125616 0.05456313  203.90356  0.000000e+00
ZZZ3         635.78304       9.312391 0.09847726   94.56387  0.000000e+00
1/2-SBSRNA4   56.97357       5.832221 0.30675574   19.01259  1.341644e-80
                     padj
                <numeric>
A1BG         1.238768e-02
A1BG-AS1     5.322729e-01
A1CF         0.000000e+00
A2LD1       2.487762e-287
A2M          0.000000e+00
...                   ...
ZYG11B       0.000000e+00
ZYX          0.000000e+00
ZZEF1        0.000000e+00
ZZZ3         0.000000e+00
1/2-SBSRNA4  1.785498e-80

 

 

ADD COMMENTlink written 19 months ago by Osvaldo20
0
gravatar for Osvaldo
19 months ago by
Osvaldo20
Spain
Osvaldo20 wrote:

opps, I've just seen your last answers.

 

Thanks, I'll download v1.11.30

ADD COMMENTlink written 19 months ago by Osvaldo20

You won't be able to run this version with the release version of R. (DESeq2 v1.12 with the fix will be released in April)

But like I said, just avoid parallel=TRUE and you will be fine.

ADD REPLYlink written 19 months ago by Michael Love14k
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 2.2.0
Traffic: 340 users visited in the last hour