Hello,

I am working with a dataset composed of three treatment groups.  I have paired samples that consist of a baseline sample and the same sample which has been enriched for active transcription.  The batch variable accounts for batches during sample preparation. The experimental design looks like this:

I need to compare C vs RT, W vs RT, and C vs W for differential expression due to treatment while controlling for the enrichment and the batch effects.   I've used design = ~Pair.nested + Batch + Treatment.  Does this design properly control for the IP/Input enrichment?  Should I approach this differently?  How should I specify the contrasts?

R version 3.1.1 (2014-07-10)

DESeq2 1.4.5

## Platform: i386-w64-mingw32/i386 (32-bit)

##

## locale:

## [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252

## [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C

## [5] LC_TIME=English_United States.1252

##

## attached base packages:

## [1] parallel stats graphics grDevices utils datasets methods base

##

## other attached packages:

## [1] RColorBrewer_1.1-2 gplots_2.16.0 DESeq2_1.4.5

## [4] RcppArmadillo_0.4.650.1.1 Rcpp_0.11.4 GenomicRanges_1.16.4

## [7] GenomeInfoDb_1.0.2 IRanges_1.22.10 BiocGenerics_0.10.0

## [10] genefilter_1.44.0 lattice_0.20-29 gtools_3.4.1

##

## loaded via a namespace (and not attached):

## [1] annotate_1.40.1 AnnotationDbi_1.26.1 Biobase_2.22.0 BiocStyle_1.0.0

## [5] bitops_1.0-6 caTools_1.17.1 DBI_0.3.1 digest_0.6.8

## [9] evaluate_0.5.5 formatR_1.0 gdata_2.13.3 geneplotter_1.40.0

## [13] grid_3.1.1 htmltools_0.2.6 KernSmooth_2.23-14 knitr_1.9

## [17] locfit_1.5-9.1 rmarkdown_0.5.1 RSQLite_1.0.0 splines_3.1.1

## [21] stats4_3.1.1 stringr_0.6.2 survival_2.37-7 tools_3.1.1

## [25] XML_3.98-1.1 xtable_1.7-4 XVector_0.2.0

hi Susan,

Looking for DE due to treatment in the enrichment samples relative to baseline can be accomplished by an interaction term. You have a 'treatment' column, but you will also need a new column, say 'type' which is Input or IP.

Then use relevel() to make sure Input is the base level (it should thanks to alphabetical order, but just to note this).

Then use a design of ~ treatment + type + treatment:pair.nested + treatment:type

Note that you can't add a batch term, because this is collinear with treatment:pair.nested. Which means that by controlling for treatment:pair.nested, you will control for differences across batch.

When running DESeq(), use betaPrior=FALSE, so that you have the standard terms (this will just make it easier for extracting contrasts with the interactions in the design).

You should then have in resultsNames(dds), [edit] two terms which look something like treatmentRT.typeIP, treatmentW.typeIP

The two terms are compared against the reference level treatment C, and so RT vs C can be extracted with:

results(dds, name="treatmentRT.typeIP")

You can compare RT vs W with 

results(dds, contrast=list("treatmentRT.typeIP","treatmentW.typeIP"))

Remember, these are *interaction terms*, so results tables are for the *extra* change in enrichment samples controlling for changes seen in Input.