I was wondering (1) how does rlog( )calculate values from raw counts, (2) the relationship between rlog counts and the effect estimator (βir), and (3) what information rlog counts tell.

In vignette:

•                     log2(qij) = Xj βi
  
  Q1: Is log2(qij) here represent the values from rlog (object)? 
  Q2: My understanding from DESeq2 paper is that log2(qij) are the values calculated from Xj βi. Is that correct?
  Q3: If Q2 is yes. What's "log qij" values for log qij=∑r Xjr βir? My understanding is: log (raw counts/size factor) instead of log2. Is that correct? 
  
  
  I saw many research made a heatmap of rlog counts. But I'm not sure what kind of information can rlog counts can tell me.
•  
• Counts comparisons from my analysis are:

• 	Control_1	Control_2	Control_3	Trt_1	Trt_2	Trt_3
  Raw count	135	258	190	0	0	0
  count(DESeq, normalized = T)	151.73	249.95	198.23	0	0	0
  rlog(object, blind=F)	6.50	6.87	6.69	5.33	5.34	5.35
  rlog(object, blind=T)	6.50	6.89	6.70	5.20	5.21	5.22
  Size factor	0.89	1.03	0.96	1.17	1.02	0.98

• And the log2fold change are:

• result(object)	lfcShrink(res, type = "normal")	lfcShrink(res, type = "apeglm")
  -10.16	-2.68	-11.24

   

  ---

  It seems either count(DESeq, normalized = T) or log2foldchange can illustrate the result better. Why we use rlog count for count heatmap? 
  
  Would you recommend to present the result by using rlog count heatmap? or log2foldchange?
  
  
  Thank you!
  Yuya

Yes, rlog is log2(q_ij). I think this formula is clear:

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#regularized-log-transformation

The formula you find in the above link is the same as in the Regularized logarithm section in the Methods of the DESeq2 paper.

rlog does not use the design matrix with covariates directly. rlog only uses the dispersion trend, and the dispersion estimates that inform the dispersion trend may (blind=FALSE) or may not (blind=TRUE) have used the design matrix.

It's perfectly fine to instead show (log or not log) normalized counts plus a pseudocount for an individual gene or for a heatmap. Normalized counts is what DESeq2 does for plotting data for individual genes for example (see plotCounts). If you pick a big enough pseudocount, log normalized counts are also good for distance calculation, our point in introducing new transformations is that the optimal pseudocount in terms of stabilizing the systematic variance trends will be different for different datasets.

What we've shown is that for calculating distances between samples over all genes the vst and rlog show some advantages over log plus a pseudocount of 1, for example, which is a common choice. So it makes sense for PCA, distance calculation, clustering. Due to the regularization, effects across samples will by definition be reduced, as they will when you add a pseudocount, but there is nevertheless an advantage for distance calculation from reducing the systematic trend in variance seen at low counts.

Note that the difference in the vst or rlog transformed counts is not the best estimator for the LFC, and this is why we do not recommend to use vst or rlog for statistical analysis of LFCs per gene in the vignette or workflow.

I ran vst and rlog on some simulated data in DESeq2:

> dds <- makeExampleDESeqDataSet(m=6)
> counts(dds)[1,] <- rep(c(0L,200L),each=3)
> counts(dds)[1,]
sample1 sample2 sample3 sample4 sample5 sample6 
      0       0       0     200     200     200 
> rld <- rlog(dds)
> assay(rld)[1,]
 sample1  sample2  sample3  sample4  sample5  sample6 
3.200444 3.215431 3.205335 6.909778 6.927411 6.884917 
> 6.9-3.2
[1] 3.7
> vsd <- varianceStabilizingTransformation(dds)
> assay(vsd)[1,]
 sample1  sample2  sample3  sample4  sample5  sample6 
4.156427 4.156427 4.156427 7.821539 7.841081 7.794033 
> 7.8 - 4.1
[1] 3.7