Hi,

I have run a DEXseq anlysis of my data, which consist of two different conditions (TG/WT) and in two different studies (S3/S5). I would like to retrieve the result tables only concerning the condition effect, where the study effect is removed. The reason why I ask for this function, is because I have worked with with the same date in a multifactor desing in DESeq2, where you provide a contrast argument. This is an example from the DESeq2 vignette:
resMFType <- results(ddsMF, contrast=c("type","single-read","paired-end"))

Is there some how an equivalent method for DEXSeq, since it is built on the same statistics and so on...

Below you can see my sample annotaions for my DEXSeqdatasat and code. I have specified two design formulae, which indicate that the StudyType factor should be treated as a blocking factor:

> sampleAnnotation(dxdC)
DataFrame with 37 rows and 5 columns
                    sample condition    study   tissue sizeFactor
                  <factor>  <factor> <factor> <factor>  <numeric>
1         18C_Ctx_z_15q_WT        WT       S3      Ctx   1.048548
2         19C_Ctx_z_15q_WT        WT       S3      Ctx   1.146222
3         20C_Ctx_z_15q_WT        WT       S3      Ctx   1.083607
4         21C_Ctx_z_15q_WT        WT       S3      Ctx   1.155945
5         22C_Ctx_z_15q_WT        WT       S3      Ctx   1.145483
...                    ...       ...      ...      ...        ...
33  148E_Ctx_Week18_15q_TG        TG       S5      Ctx  0.9056882
34  149E_Ctx_Week18_15q_TG        TG       S5      Ctx  0.9465755
35  150E_Ctx_Week18_15q_TG        TG       S5      Ctx  0.7392238
36  151E_Ctx_Week18_15q_TG        TG       S5      Ctx  0.9036385
37  152E_Ctx_Week18_15q_TG        TG       S5      Ctx  1.1340215

formulaFullModel    =  ~ sample + exon + study:exon + condition:exon
formulaReducedModel =  ~ sample + exon + study:exon

dxdC = estimateDispersions( dxdC, BPPARAM=BPPARAM, formula = formulaFullModel )

dxdC = testForDEU( dxdC,
        reducedModel = formulaReducedModel,
        fullModel = formulaFullModel, BPPARAM=BPPARAM )
        
dxrC = DEXSeqResults( dxdC )
 

Thanks a lot for your time,

Best, Maria Dalby

> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] C

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

other attached packages:
 [1] DEXSeq_1.12.1             BiocParallel_1.0.1       
 [3] DESeq2_1.6.3              RcppArmadillo_0.4.600.4.0
 [5] Rcpp_0.11.4               GenomicRanges_1.18.4     
 [7] GenomeInfoDb_1.2.4        IRanges_2.0.1            
 [9] S4Vectors_0.4.0           Biobase_2.26.0           
[11] BiocGenerics_0.12.1      

loaded via a namespace (and not attached):
 [1] AnnotationDbi_1.28.1 BBmisc_1.8           BatchJobs_1.5       
 [4] Biostrings_2.34.1    DBI_0.3.1            Formula_1.2-0       
 [7] Hmisc_3.14-6         MASS_7.3-37          RColorBrewer_1.1-2  
[10] RCurl_1.95-4.5       RSQLite_1.0.0        Rsamtools_1.18.2    
[13] XML_3.98-1.1         XVector_0.6.0        acepack_1.3-3.3     
[16] annotate_1.44.0      base64enc_0.1-2      biomaRt_2.22.0      
[19] bitops_1.0-6         brew_1.0-6           checkmate_1.5.1     
[22] cluster_1.15.3       codetools_0.2-10     colorspace_1.2-4    
[25] digest_0.6.8         fail_1.2             foreach_1.4.2       
[28] foreign_0.8-62       genefilter_1.48.1    geneplotter_1.44.0  
[31] ggplot2_1.0.0        grid_3.1.2           gtable_0.1.2        
[34] hwriter_1.3.2        iterators_1.0.7      lattice_0.20-29     
[37] latticeExtra_0.6-26  locfit_1.5-9.1       munsell_0.4.2       
[40] nnet_7.3-8           plyr_1.8.1           proto_0.3-10        
[43] reshape2_1.4.1       rpart_4.1-8          scales_0.2.4        
[46] sendmailR_1.2-1      splines_3.1.2        statmod_1.4.20      
[49] stringr_0.6.2        survival_2.37-7      tools_3.1.2         
[52] xtable_1.7-4         zlibbioc_1.12.0     

The code you posted is correct. It will give you all genes whose expression differs by condition, after blocking for study. Is this not what you want?

If, just out of curiosity, you want to know whether there is a differences between studies, you can, of course, run a second analysis, this time using condition as blocking factor ans study as factor of interest.

The reason that you have to run two separate analyses in DEXSeq, but can get p values for both question from the same run in DESeq2, is that in DEXSeq we always use a likelihood ratio test (LRT), which compares two models, which should differ in only one factor. In DESeq2, we fit (by default) only the full model and then perform Wald tests, which can return a p value for each coefficient in the model, without the need to fit a reduced model.

(There is lot to be said about the pros and cons of LRT vs Wald test in this context, but this would become a rather lengthy post.)