Search
Question: Question: DESeq2 - factor base level changes LFCs -- hoping to understand this better
0
gravatar for Victor Missirian
21 months ago by
USA/Davis/UC Davis
Victor Missirian0 wrote:

Hi Michael,

    I'd like to ask a question related to a previous post: "Question: DESeq2 - factor base level changes the DE genes."

    I ran the following code on DESeq2 (version 1.4.5):

## START CODE

sampleTable <- read.delim(sample.table.file, row.names=1, stringsAsFactors=FALSE)
coldata <- DataFrame(sample=sampleTable$SampleName,
            group=factor(sampleTable$group),
            condition=factor(sampleTable$condition))

countdata <- read.delim(htseq.output.file, header=TRUE, sep="\t", quote="", comment.char="", row.names=1)


## using "ref.group" as reference level for "group" factor

dds.1 <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ group + condition + group:condition)
dds.1$condition <- relevel(dds.1$condition, ref.condition)
dds.1$group <- relevel(dds.1$group, ref.group)
dds.1 <- DESeq(dds.1)
res.other.vs.ref.condition.effect.in.ref.group <- results(dds.1, contrast=c("condition", other.condition, ref.condition))
res.ratio.of.other.vs.ref.condition.effect.in.other.vs.ref.group <- results(dds.1, name=paste("group", other.group, ".condition", other.condition, sep=""))


## using "other.group" as reference level for "group" factor
dds.2 <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ group + condition + group:condition)
dds.2$condition <- relevel(dds.2$condition, ref.condition)
dds.2$group <- relevel(dds.2$group, other.group)
dds.2 <- DESeq(dds.2)
res.other.vs.ref.condition.effect.in.other.group <- results(dds.2, contrast=c("condition", other.condition, ref.condition))

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] hash_3.0.1                DESeq2_1.4.5              RcppArmadillo_0.5.100.1.0
 [4] Rcpp_0.11.6               GenomicRanges_1.16.4      GenomeInfoDb_1.0.2       
 [7] IRanges_1.22.10           Biobase_2.24.0            BiocGenerics_0.10.0      
[10] pasilla_0.4.0            

loaded via a namespace (and not attached):
 [1] annotate_1.42.1      AnnotationDbi_1.26.1 DBI_0.3.1            DESeq_1.16.0        
 [5] genefilter_1.46.1    geneplotter_1.42.0   grid_3.1.0           lattice_0.20-31     
 [9] locfit_1.5-9.1       RColorBrewer_1.1-2   RSQLite_1.0.0        splines_3.1.0       
[13] stats4_3.1.0         survival_2.38-1      tools_3.1.0          XML_3.98-1.1        
[17] xtable_1.7-4         XVector_0.4.0  

## END CODE

 

And I found that the LFC of the interaction effect was not the same as the difference between the LFCs of the condition effects for the two groups:

> res.other.vs.ref.condition.effect.in.ref.group["AT5G20230","log2FoldChange"]
[1] 1.632546
> res.ratio.of.other.vs.ref.condition.effect.in.other.vs.ref.group["AT5G20230","log2FoldChange"]
[1] 2.30044
> res.other.vs.ref.condition.effect.in.other.group["AT5G20230","log2FoldChange"]
[1] 4.312705


    If possible, could you please explain the reason that this is the case?

    Do the values computed for the mean normalized read counts or dispersions depend on the reference level of the "group" factor?

    If so, is there a quick explanation that could help me to better understand how this dependence works?


Thanks,
-Victor Missirian

ADD COMMENTlink modified 21 months ago by Michael Love12k • written 21 months ago by Victor Missirian0
0
gravatar for Michael Love
21 months ago by
Michael Love12k
United States
Michael Love12k wrote:

The reason is that for DESeq2 versions 1.4-1.8 the main effect terms are not shrunk, while the interaction terms are shrunk. So in comparison 1 you are adding (main + interaction) and you have (unshrunk + shrunk) while in comparison 2 you have just a main effect terms which is not shrunk.

I recommend turning off the betaPrior for designs with interactions, and in the next release (v1.10 in Oct 2015), it will automatically turn off the betaPrior. This looks like:

dds <- DESeq(dds, betaPrior=FALSE)

For more background, I wrote this other post recently on the topic:

A: Interactions in DESeq2

ADD COMMENTlink modified 21 months ago • written 21 months ago by Michael Love12k

Thanks!  I will try turning off betaPrior!

-Victor

ADD REPLYlink written 21 months ago by Victor Missirian0

In R 1.8.1 I have done a small test for a design with no interactions of the type ~line+treatment, where line has 4 different levels.

I have noticed (small) differences in the estimated pvalues and adjusted pvalues if I change names of the levels for "line" factor (so that the first level in alphabetical order that appears to be used as reference is going to change).

This is true for betaPRIOR=FALSE as well.

I was curious to know whether these differences are because of how estimation of parameters is implemented? And if yes, are they expected to be always minor?

Many thanks

ADD REPLYlink modified 12 months ago • written 12 months ago by deut1340

What version of DESeq2 are you using?

ADD REPLYlink written 12 months ago by Michael Love12k

Sorry, I meant to write DESeq2 1.8.1 - I have to say I have not (yet) tried with the current 1.12.3

ADD REPLYlink written 12 months ago by deut1340

You should only ever see very small differences (not in the first significant digits) having to do with algorithm convergence.

ADD REPLYlink written 12 months ago by Michael Love12k

Thanks for the clarification!

ADD REPLYlink written 12 months ago by deut1340
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: 114 users visited in the last hour