[BioC] DESeq2 - factor base level changes the DE genes.

[Ximena [guest]]

Hello,

I have an RNAseq experiment from two strains of mice, that have been exposed to two different environments (three biological replicates of each combination). So my experimental setup is

> colData
        genetics environment
black2        B6       black
black3        B6       black
black5        B6       black
black1        B6      yellow
black4        B6      yellow
black6        B6      yellow
agouti1      129      yellow
agouti4      129      yellow
agouti6      129      yellow
agouti2      129       black
agouti3      129       black
agouti5      129       black

I want to test which genes change their expression when the environment changes, controlling for the genetic background. I am using DESeq2 for this, with the formula ~ genetics + environment + genetics:environment.

However, the choice of the base level for the factors in my experimental design changes the genes that are significant. If I leave the default reference levels (129 for genetics and black for env) I get more DE genes that if I set as reference B6 for genetics and yellow for environemnt.
I thought the reference level is used to decide how to calculate the fold change and is necessary to correctly interpret the coefficients, but shouldn't the results be the same?

Any thoughts on what am I missing?

Thank you very much in advance.


My code is below.

## 1 ##
> levels(colData$genetics)
[1] "129" "B6" 
> levels(colData$environment)
[1] "black"  "yellow"

> dds <- DESeqDataSetFromMatrix(countData = data[,8:19], colData = colData, design = ~ genetics + environment + genetics:environment)
> dds <- DESeq(dds)

> resG <- results(dds, alpha=0.05, name="genetics_B6_vs_129")
> resE <- results(dds, alpha=0.05, name="environment_yellow_vs_black")
> resGE <- results(dds, alpha=0.05, name="geneticsB6.environmentyellow")

> dim(subset(resG, resG$padj < 0.05))
[1] 3675    6
> dim(subset(resE, resE$padj < 0.05))
[1] 20  6
> dim(subset(resGE, resGE$padj < 0.05))
[1] 24  6

## 2 ##
> colData2 <- colData
> colData2$genetics <- relevel(colData2$genetics, ref="B6")
> colData2$environment <- relevel(colData2$environment, ref="yellow")
> levels(colData2$genetics)
[1] "B6"  "129"
> levels(colData2$environment)
[1] "yellow" "black" 

> dds2 <- DESeqDataSetFromMatrix(countData = data[,8:19], colData = colData2, design = ~ genetics + environment + genetics:environment)
> dds2 <- DESeq(dds2)

> res2G <- results(dds2, alpha=0.05, name="genetics_129_vs_B6")
> res2E <- results(dds2, alpha=0.05, name="environment_black_vs_yellow")
> res2GE <- results(dds2, alpha=0.05, name="genetics129.environmentblack")

> dim(subset(res2G, resG$padj < 0.05))
[1] 3127    6
> dim(subset(res2E, resE$padj < 0.05))
[1] 2 6
> dim(subset(res2GE, resGE$padj < 0.05))
[1] 0 6

##

> g1 <- subset(resG, resG$padj < 0.05)
> g2 <- subset(res2G, res2G$padj < 0.05)

> length(intersect(rownames(g1), rownames(g2)))
[1] 2640



 -- output of sessionInfo(): 

R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

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

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

other attached packages:
[1] DESeq2_1.2.10           RcppArmadillo_0.4.320.0 Rcpp_0.11.2            
[4] GenomicRanges_1.14.4    XVector_0.2.0           IRanges_1.20.7         
[7] BiocGenerics_0.8.0     

loaded via a namespace (and not attached):
 [1] annotate_1.40.1      AnnotationDbi_1.24.0 Biobase_2.22.0       DBI_0.2-7           
 [5] genefilter_1.44.0    grid_3.0.2           lattice_0.20-23      locfit_1.5-9.1      
 [9] RColorBrewer_1.0-5   RSQLite_0.11.4       splines_3.0.2        stats4_3.0.2        
[13] survival_2.37-7      tools_3.0.2          XML_3.95-0.2         xtable_1.7-3 

--
Sent via the guest posting facility at bioconductor.org.