Answers.Body: <p>It depends what you mean by “normalised”. As Devon said, the <code>normalized = TRUE</code> argument to the <code>count</code> function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the <code>sizeFactors(dds)</code>).</p>
  
  <p>However, <a href="https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization" rel="noreferrer">as the vignette explains</a>, downstream processing generally requires more advanced normalisation, to account for the <a href="https://en.wikipedia.org/wiki/Heteroscedasticity" rel="noreferrer">heteroscedasticity</a> of the counts. This is often done by simply <code>log</code>ging the counts but this has obvious drawbacks (most trivially, what do we do with 0 counts? A workaround is to add a pseudocount but that’s problematic too).</p>
  
  <p>DESeq2 offers two different methods to perform a more rigorous analysis:</p>
  
  <ul>
  <li><code>rlog</code> — a <em>regularised log</em>, and</li>
  <li><code>vst</code> — a <em>variance stabilising transformation</em>.</li>
  </ul>
  
  <p>You’d generally use either of these for downstream analysis, not <code>count(dds, normalized = TRUE)</code>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-29T05:44:28.210Z
Answers.Text: to make my confusion clear let's say im going for any parametric or non parametric test between any comparison for genes ,can i take the rlog or vst ? As rlog and vst i take it for clustering ,pca  ,correlation ,wgcna .So is it possible ?
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-29T12:13:52.663Z
Answers.Text: @krushnachChandra It depends on the test and the comparison. DESeq for instance does *not* use either rlog or vst for its test of differential expression, as these would be inappropriate.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-29T12:25:53.243Z
Answers.Text: yes that i read in the paper , for example i got a cluster and i wanted to see what if the genes are having significant difference in their expression so if i have to run say a non parametric test can i use those rlog ? values or i have to use the count normalised ones ?
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-29T13:46:15.223Z
Answers.Text: @krushnachChandra You cannot use count normalised values in this scenario, nor simple rlog values, since they are not normalised for transcript length (and thus longer transcripts will show higher expression on average). You could use something like size-factor normalised [TPM](https://bioinformatics.stackexchange.com/a/69/29) instead. You can also log those TPMs but for a nonparametric test it would make no difference — that’s kind of the point of non-parametric tests.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-30T06:10:50.640Z
Answers.Text: So for PCA ,clustering i can use rlog i hope that part im doing correct let me know, i read your post , so can i use the deseq2 value and convert them into TPM ,but there is one part where even deseq2 also calculate size factor "dds <- estimateSizeFactors(dds)" this i guess . So can i use your TPM conversion method on deseq2 values? or just counts .
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-30T11:07:59.420Z
Answers.Text: im curious about this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4986529/ Figure 6 last box plot "(Right) Mean expression of the two gene sets in each population. Error bars represent SD. CPM, counts per million" where they use cpm values i  guess these are edge R output ,so here they used that value to make comparison ,so i was wondering if it can be done with rlog values ?sorry for mundane question my concept about these are very fuzzy as im not a statistician ..
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-24T10:17:43.520Z
Answers.LastActivityDate: 2017-05-24T10:17:43.520Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>The normalized counts themselves can be accessed with <code>counts(dds, normalized=T)</code>.</p>
  
  <p>Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that <code>geno_treat</code> has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in <code>geno_treat</code>.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-24T07:44:58.903Z
Answers.Text: `mean(assays(dds)$mu["WBGene00000001",] / sizeFactors(dds))` gives something very close to `res["WBGene00000001",]$baseMean` (181.9139). The vignette says that `assays(dds)$mu` correspond to s_j * q_ij where, according to the DESeq2 paper, s_j is the size factor for library j and q_ij is, if I understand correctly, an estimation of a quantity "proportional to the concentration of cDNA fragments from the gene in the sample". Is this q_ij what you call "base factor" ?
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-24T09:52:59.093Z
Answers.Text: Not really. With an "expanded model matrix" the intercept (aka, "baseMean") should be approximately the average of the normalized counts. For non-"expanded model matrices" then one level of the `group` factor (the "base level") becomes the intercept, so the baseMean is then just the average of its counts. The base level is the alphanumerically first level in R by default, if you have `grp = factor(c("WT", "MUT"))` then the `MUT` group is used to calculate the baseMean (assuming no expanded model matrices).
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-24T12:13:44.303Z
Answers.Text: Here, I find that the mean of all normalized values (`mean(counts(dds, normalized=T)["WBGene00000001",])`), is actually equal to the baseMean (`res["WBGene00000001",]$baseMean`). This is not the case if I compute the mean only for a given level of the `geno_treat` factor.
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-24T12:58:41.690Z
Answers.Text: I guess it using an "expanded model matrix" then. I know Michael used to code in exactly the cases when that would happen, but perhaps he's since generalized the code, which is nice :)
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-22T18:17:38.167Z
Answers.LastActivityDate: 2017-05-22T18:17:38.167Z
Answers.UserDisplayName: Devon Ryan
Body: <p>The results obtained by running the <code>results</code> command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.</p>
  
  <p><strong>How can I access the normalized counts proper?</strong></p>
  
  <p>I tried the following (continuing with the example used <a href="https://bioinformatics.stackexchange.com/a/192/292">here</a>):</p>
  
  <pre><code>&gt; dds &lt;- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
  &gt; dds &lt;- DESeq(dds)
  estimating size factors
  estimating dispersions
  gene-wise dispersion estimates
  mean-dispersion relationship
  final dispersion estimates
  fitting model and testing
  &gt; res &lt;- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))
  </code></pre>
  
  <p>Here is what I have for the first gene:</p>
  
  <pre><code>&gt; res["WBGene00000001",]$baseMean
  [1] 181.7862
  &gt; mean(assays(dds)$mu["WBGene00000001",])
  [1] 231.4634
  &gt; mean(assays(dds)$counts["WBGene00000001",])
  [1] 232.0833
  </code></pre>
  
  <p><code>assays(dds)$counts</code> corresponds to the raw counts. <code>assays(dds)$mu</code> seems to be a transformation of these counts approximately preserving their mean, but this mean is very different from the "baseMean" value, so these are likely not the normalized values.</p>
  
Tags: <r><bioconductor><deseq2><normalization>
Title: How can I extract normalized read count values from DESeq2 results?
ThreadText: # How can I extract normalized read count values from DESeq2 results?
  
  (Asked by: username_0 on 2017-05-22)
  
  The results obtained by running the `results` command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.
  
  **How can I access the normalized counts proper?**
  
  I tried the following (continuing with the example used [here](https://bioinformatics.stackexchange.com/a/192/292)):
  
  ```
  > dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
  > dds <- DESeq(dds)
  estimating size factors
  estimating dispersions
  gene-wise dispersion estimates
  mean-dispersion relationship
  final dispersion estimates
  fitting model and testing
  > res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))
  
  ```
  
  Here is what I have for the first gene:
  
  ```
  > res["WBGene00000001",]$baseMean
  [1] 181.7862
  > mean(assays(dds)$mu["WBGene00000001",])
  [1] 231.4634
  > mean(assays(dds)$counts["WBGene00000001",])
  [1] 232.0833
  
  ```
  
  `assays(dds)$counts` corresponds to the raw counts. `assays(dds)$mu` seems to be a transformation of these counts approximately preserving their mean, but this mean is very different from the "baseMean" value, so these are likely not the normalized values.
  
  **Tags:** r, bioconductor, deseq2, normalization
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2017-05-22. Score: 13 | Accepted answer
  
  The normalized counts themselves can be accessed with `counts(dds, normalized=T)`.
  
  Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that `geno_treat` has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in `geno_treat`.
  
  > **Comments:**
  > **username_0**: `mean(assays(dds)$mu["WBGene00000001",] / sizeFactors(dds))` gives something very close to `res["WBGene00000001",]$baseMean` (181.9139). The vignette says that `assays(dds)$mu` correspond to s_j * q_ij where, according to the DESeq2 paper, s_j is the size factor for library j and q_ij is, if I understand correctly, an estimation of a quantity "proportional to the concentration of cDNA fragments from the gene in the sample". Is this q_ij what you call "base factor" ?  
  > **username_3**: Not really. With an "expanded model matrix" the intercept (aka, "baseMean") should be approximately the average of the normalized counts. For non-"expanded model matrices" then one level of the `group` factor (the "base level") becomes the intercept, so the baseMean is then just the average of its counts. The base level is the alphanumerically first level in R by default, if you have `grp = factor(c("WT", "MUT"))` then the `MUT` group is used to calculate the baseMean (assuming no expanded model matrices).  
  > **username_0**: Here, I find that the mean of all normalized values (`mean(counts(dds, normalized=T)["WBGene00000001",])`), is actually equal to the baseMean (`res["WBGene00000001",]$baseMean`). This is not the case if I compute the mean only for a given level of the `geno_treat` factor.  
  > **username_3**: I guess it using an "expanded model matrix" then. I know Michael used to code in exactly the cases when that would happen, but perhaps he's since generalized the code, which is nice :)  
  
  ---
  
  ### Answer by username_1 on 2017-05-24. Score: 6
  
  It depends what you mean by “normalised”. As username_3 said, the `normalized = TRUE` argument to the `count` function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the `sizeFactors(dds)`).
  
  However, [as the vignette explains](https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization), downstream processing generally requires more advanced normalisation, to account for the [heteroscedasticity](https://en.wikipedia.org/wiki/Heteroscedasticity) of the counts. This is often done by simply `log`ging the counts but this has obvious drawbacks (most trivially, what do we do with 0 counts? A workaround is to add a pseudocount but that’s problematic too).
  
  DESeq2 offers two different methods to perform a more rigorous analysis:
  
  * `rlog` — a *regularised log*, and
  * `vst` — a *variance stabilising transformation*.
  
  You’d generally use either of these for downstream analysis, not `count(dds, normalized = TRUE)`.
  
  > **Comments:**
  > **username_2**: to make my confusion clear let's say im going for any parametric or non parametric test between any comparison for genes ,can i take the rlog or vst ? As rlog and vst i take it for clustering ,pca  ,correlation ,wgcna .So is it possible ?  
  > **username_1**: @krushnachChandra It depends on the test and the comparison. DESeq for instance does *not* use either rlog or vst for its test of differential expression, as these would be inappropriate.  
  > **username_2**: yes that i read in the paper , for example i got a cluster and i wanted to see what if the genes are having significant difference in their expression so if i have to run say a non parametric test can i use those rlog ? values or i have to use the count normalised ones ?  
  > **username_1**: @krushnachChandra You cannot use count normalised values in this scenario, nor simple rlog values, since they are not normalised for transcript length (and thus longer transcripts will show higher expression on average). You could use something like size-factor normalised [TPM](https://bioinformatics.stackexchange.com/a/69/29) instead. You can also log those TPMs but for a nonparametric test it would make no difference — that’s kind of the point of non-parametric tests.  
  > **username_2**: So for PCA ,clustering i can use rlog i hope that part im doing correct let me know, i read your post , so can i use the deseq2 value and convert them into TPM ,but there is one part where even deseq2 also calculate size factor "dds <- estimateSizeFactors(dds)" this i guess . So can i use your TPM conversion method on deseq2 values? or just counts .  
  > **username_2**: im curious about this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4986529/ Figure 6 last box plot "(Right) Mean expression of the two gene sets in each population. Error bars represent SD. CPM, counts per million" where they use cpm values i  guess these are edge R output ,so here they used that value to make comparison ,so i was wondering if it can be done with rlog values ?sorry for mundane question my concept about these are very fuzzy as im not a statistician ..  
  
  ---
  
  