DESEQ2: apparently wrong log2(fold-change) calculation
5
0
Entering edit mode
Osvaldo ▴ 30
@osvaldo-6150
Last seen 13 months ago
Spain

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 !

 

 

 

deseq2 differential gene expression log2fc • 4.4k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 1 hour ago
United States
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 COMMENT
0
Entering edit mode
Osvaldo ▴ 30
@osvaldo-6150
Last seen 13 months ago
Spain

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 COMMENT
0
Entering edit mode
What is printed above res when you print to console?
ADD REPLY
0
Entering edit mode
Osvaldo ▴ 30
@osvaldo-6150
Last seen 13 months ago
Spain

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

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

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

Fixed in v1.11.30

ADD REPLY
0
Entering edit mode
Osvaldo ▴ 30
@osvaldo-6150
Last seen 13 months ago
Spain

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 COMMENT
0
Entering edit mode
Osvaldo ▴ 30
@osvaldo-6150
Last seen 13 months ago
Spain

opps, I've just seen your last answers.

 

Thanks, I'll download v1.11.30

ADD COMMENT
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 531 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