Clarify design formula and blind dispersion estimation

Author: pmoris

bin/deseq2_qc.r

@@ -92,7 +92,8 @@ if (decompose) {
 DDSFile <- paste(opt$outprefix,".dds.RData",sep="")
 
 counts  <- count.table[,samples.vec,drop=FALSE]
-dds     <- DESeqDataSetFromMatrix(countData=round(counts), colData=coldata, design=~ 1)
+# `design=~1` creates intercept-only model, equivalent to setting `blind=TRUE` for transformation.
+dds     <- DESeqDataSetFromMatrix(countData=round(counts), colData=coldata, design=~1)
 dds     <- estimateSizeFactors(dds)
 if (min(dim(count.table))<=1)  { # No point if only one sample, or one gene
     save(dds,file=DDSFile)
@@ -102,10 +103,10 @@ if (min(dim(count.table))<=1)  { # No point if only one sample, or one gene
 }
 if (!opt$vst) {
     vst_name <- "rlog"
-    rld      <- rlog(dds)
+    rld      <- rlog(dds, blind=TRUE)   # blind=TRUE is the default and already implied by design=~1
 } else {
     vst_name <- "vst"
-    rld      <- varianceStabilizingTransformation(dds)
+    rld      <- varianceStabilizingTransformation(dds, blind=TRUE)
 }
 
 assay(dds, vst_name) <- assay(rld)

docs/output.md

@@ -642,6 +642,8 @@ The script included in the pipeline uses DESeq2 to normalise read counts across
 
 By default, the pipeline uses the `vst` transformation which is more suited to larger experiments. You can set the parameter `--deseq2_vst false` if you wish to use the DESeq2 native `rlog` option. See [DESeq2 docs](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization) for a more detailed explanation.
 
+Both types of transformation are performed blind, i.e. using across-all-samples variability, without using any prior information on experimental groups (equivalent to using an intercept-only design), as recommended by the [DESeq2 docs](https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#blind-dispersion-estimation).
+
 The PCA plots are generated based alternately on the top five hundred most variable genes, or all genes. The former is the conventional approach that is more likely to pick up strong effects (ie the biological signal) and the latter, when different, is picking up a weaker but consistent effect that is synchronised across many transcripts. We project both of these onto the first two PCs (shown in the top row of the figure below), which is the best two dimensional representation of the variation between samples.
 
 We also explore higher components in terms of experimental factors inferred from sample names. If your sample naming convention follows a strict policy of using underscores to delimit values of experimental factors (for example `WT_UNTREATED_REP1`) and all names have the same number of underscores (so excluding `WT_TREATED_10ml_REP1` from being compatible with the previous label), then any of these factors that are informative (ie label some but not all samples the same) then we individually plot upto the first five PCs, per experimental level, for each of the experimental factors.