hi Richard,

On Mon, May 27, 2013 at 5:11 PM, Hickman, R.J. (Richard)
<R.J.Hickman@uu.nl>wrote:

>
>
> 2) For genes that are strongly expressed, the rLog transformation does not
> perform how I would expect, in that it seems to boost the values of the
> untreated samples and reduce the values of the treated samples (see
> attached figure "DESeq_transformations_timeseries_strongExpression.png"),
> meaning that the fold changes resulting from the transformed data are
> greatly reduced. As these expression levels are what I would class as
> "high", then I would have thought that these values should remain
> relatively unchanged by the transformation.
>
>
The rlogTransformation is squashing the values here because it does not
know that there are biological replicates (and so the dispersion should be
small), it just sees altogether the 2 series * 10 time points * 4
replicates which show a lot of variance. I guess then in this case I would
not recommend unsupervised dispersion estimation.

Two approaches to get a transformation which "sees" the replicate structure
are: rlogTransformation() has an argument 'samplesVector', which allows you
to fit a single value for each sample if there are replicates. However I
think you want to see the variance at each time point, so another approach
would be to obtain a dispersion estimates which "sees" the replicates.

With some simulated data, I try to show the problem you see and a solution.
 Here I simulate 5 time points with 4 replicates each:

> library(DESeq2)
> options(digits=2)
> dds <- makeExampleDESeqDataSet(m=20)
> colData(dds)$time <- factor(rep(1:5,each=4))
> colnames(dds) <- paste0("t",rep(1:5,each=4),"-",1:4)

Make the first row have increasing mean values: 2^5, 2^6, 2^7, 2^8, 2^9

> counts(dds)[1,] <- rpois(20, 2^rep(c(5:9),each=4))
> counts(dds)[1,]
t1-1 t1-2 t1-3 t1-4 t2-1 t2-2 t2-3 t2-4 t3-1 t3-2 t3-3 t3-4 t4-1 t4-2 t4-3
t4-4 t5-1 t5-2 t5-3
  34   30   35   28   72   61   57   61  130  148  123  116  245  250  285
 253  506  552  454
t5-4
 507

With DESeq2 1.0.10, we added an argument 'unsupervised', with default TRUE,
which internally replaces design(dds) <- ~ 1 and re-estimates dispersions.
 Here the values are squashed together because the row has a high
dispersion estimate.

> rld <- rlogTransformation(dds, unsupervised=TRUE)
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
> assay(rld)[1,]
t1-1 t1-2 t1-3 t1-4 t2-1 t2-2 t2-3 t2-4 t3-1 t3-2 t3-3 t3-4 t4-1 t4-2 t4-3
t4-4 t5-1 t5-2 t5-3
 6.4  6.3  6.4  6.3  6.7  6.6  6.6  6.6  7.0  7.1  7.0  7.0  7.5  7.5  7.6
 7.5  8.1  8.2  8.0
t5-4
 8.1

But if you want the transformation to "see" the replicate structure, you
can supply the time points like so and estimate dispersions:

> design(dds) <- ~ time
> dds <- estimateSizeFactors(dds)
> dds <- estimateDispersions(dds)
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

Then specify 'unsupervised' to be FALSE:

> rld <- rlogTransformation(dds, unsupervised=FALSE)

Now the transformed values for this row are nearly the same as log2 values:

> assay(rld)[1,]
t1-1 t1-2 t1-3 t1-4 t2-1 t2-2 t2-3 t2-4 t3-1 t3-2 t3-3 t3-4 t4-1 t4-2 t4-3
t4-4 t5-1 t5-2 t5-3
 5.2  5.0  5.2  5.0  6.1  6.0  5.9  5.9  7.0  7.2  6.9  6.8  7.8  7.8  8.0
 7.9  8.8  8.9  8.6
t5-4
 8.8
> log2(counts(dds[1,]))
         t1-1 t1-2 t1-3 t1-4 t2-1 t2-2 t2-3 t2-4 t3-1 t3-2 t3-3 t3-4 t4-1
t4-2 t4-3 t4-4 t5-1
feature1  5.1  4.9  5.1  4.8  6.2  5.9  5.8  5.9    7  7.2  6.9  6.9  7.9
 8  8.2    8    9
         t5-2 t5-3 t5-4
feature1  9.1  8.8    9


best,

Mike

	[[alternative HTML version deleted]]

