Hi,

Can someone explain why the number of differentially expressed genes at FDR<10% changes when the order of variables changes in the design?

Option1: Condition given at the end of the design matrix -  54 genes at FDR<10%

>dds=DESeqDataSetFromMatrix(countData=countData,colData=coldata,design=~Race+RIN_Acsg+GCpercent+condition)

> dds=DESeq(dds)                                    

> ctrlBP=results(dds,contrast=c("condition","Bipolar","Control"))   

> ctrlBP=ctrlBP[order(ctrlBP$padj),]                                
> summary(ctrlBP)                                                   
                                                                    
out of 21228 with nonzero total read count                          
adjusted p-value < 0.1                                              
LFC > 0 (up)     : 30, 0.14%                                        
LFC < 0 (down)   : 24, 0.11%                                        
outliers [1]     : 0, 0%                                            
low counts [2]   : 0, 0%                                            
(mean count < 0)                                                    
[1] see 'cooksCutoff' argument of ?results                          
[2] see 'independentFiltering' argument of ?results                 

Option2: Condition given at the beginning of the design matrix -  24 genes at FDR<10%

>dds=DESeqDataSetFromMatrix(countData=countData,colData=coldata,design=~condition+Race+RIN_Acsg+GCpercent)

> dds=DESeq(dds)                                     

> ctrlBP=results(dds,contrast=c("condition","Bipolar","Control"))   

> ctrlBP=ctrlBP[order(ctrlBP$padj),]                    

> summary(ctrlBP)                                                 

out of 21228 with nonzero total read count

adjusted p-value < 0.1

LFC > 0 (up)       : 13, 0.061%

LFC < 0 (down)     : 11, 0.052%

outliers [1]       : 0, 0%

low counts [2]     : 0, 0%

(mean count < 0)

[1] see 'cooksCutoff' argument of ?results

[2] see 'independentFiltering' argument of ?results

Thanks,

Nirmala

Can you run this again to double check? This example has identical p-values:

> dds <- makeExampleDESeqDataSet()
> dds$batch <- factor(rep(c(1,2,1,2),each=3))
> design(dds) <- ~batch + condition
> dds <- DESeq(dds)
> res <- results(dds, name="condition_B_vs_A")
> design(dds) <- ~condition + batch
> dds <- DESeq(dds)
> res2 <- results(dds, name="condition_B_vs_A")
> all.equal(res$pvalue, res2$pvalue)
[1] TRUE