Answers.Body: <p>You could correct your logCPMs (<code>edgeR::cpm(y, log=TRUE)</code>) with the covariates from RUVseq using limma and then use these for plotting:</p>
  <pre><code>set2 &lt;- RUVg(...)
  
  designRUV &lt;- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))
  
  y &lt;- DGEList(counts=counts(set2), group=pheno$condition4)
  y &lt;- calcNormFactors(y, method=&quot;upperquartile&quot;)
  y &lt;- estimateGLMCommonDisp(y, designRUV)
  y &lt;- estimateGLMTagwiseDisp(y, designRUV)
  fit &lt;- glmFit(y, designRUV)
  
  # correction with limma
  logcpm &lt;- cpm(y, log=TRUE)
  design &lt;- model.matrix(~0+condition4, pData(set2))
  corrected &lt;- removeBatchEffect(logcpm, design=design, covariates=as.matrix(pData(set2)[,c(W_1,&quot;W_2&quot;])
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-07-19T06:46:05.523Z
Answers.Text: Thank you @alderaan_shot_first. Can I ask why you used a "-" in the `model.matrix`? Also, after reading your answer I have been checking the `removeBatchEffect` documentation and when you put the covariates argument you are setting all the phenoData of set2... Don't would be `covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])` ?. Thank you!!!
Answers.UserDisplayName: FrAoJm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-07-19T12:00:08.643Z
Answers.Text: was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit.
Answers.UserDisplayName: alderaan_shot_first
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-07-18T22:34:44.660Z
Answers.LastActivityDate: 2022-07-19T18:19:41.650Z
Answers.LastEditDate: 2022-07-19T18:19:41.650Z
Answers.UserDisplayName: alderaan_shot_first
Body: <p>I have used the empirical RUVg method (from RUVseq) to estimate the unwanted variation of my dataset (consisting of several public datasets analysed together, with controls and case samples in different datasets.</p>
  <p>Afterwards, I used EdgeR to perform my DE analysis with 2 RUV covariates:</p>
  <pre><code>designRUV &lt;- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))
  
  y &lt;- DGEList(counts=counts(set2), group=pheno$condition4)
  y &lt;- calcNormFactors(y, method=&quot;upperquartile&quot;)
  y &lt;- estimateGLMCommonDisp(y, designRUV)
  y &lt;- estimateGLMTagwiseDisp(y, designRUV)
  fit &lt;- glmFit(y, designRUV)
  </code></pre>
  <p>I got my results after making all the contrast that I needed but I would like to plot the expression of some genes of interest. One plot per gene with all the conditions in <code>condition4</code> has 9 categories (expression level or norm counts on the Y axis).</p>
  <p>My doubt is which counts to use. I have the normalised ones from RUVg, but it states they should only use them for exploratory analysis (maybe this is the case??). These counts as I understand are the residuals of the model, so not really sure what to do.</p>
  <p>Another option is to use the normalised counts produced by EdgeR, but these ones doesn't account for the unwanted variation (W_1, W_2 covariates) (right???)</p>
  <p>Any advice or nice solution?? :=)
  Thanks in advance!!</p>
  
Tags: <rna-seq><normalization><edger>
Title: plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
ThreadText: # plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
  
  (Asked by: username_0 on 2022-07-18)
  
  I have used the empirical RUVg method (from RUVseq) to estimate the unwanted variation of my dataset (consisting of several public datasets analysed together, with controls and case samples in different datasets.
  
  Afterwards, I used EdgeR to perform my DE analysis with 2 RUV covariates:
  
  ```
  designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))
  
  y <- DGEList(counts=counts(set2), group=pheno$condition4)
  y <- calcNormFactors(y, method="upperquartile")
  y <- estimateGLMCommonDisp(y, designRUV)
  y <- estimateGLMTagwiseDisp(y, designRUV)
  fit <- glmFit(y, designRUV)
  
  ```
  
  I got my results after making all the contrast that I needed but I would like to plot the expression of some genes of interest. One plot per gene with all the conditions in `condition4` has 9 categories (expression level or norm counts on the Y axis).
  
  My doubt is which counts to use. I have the normalised ones from RUVg, but it states they should only use them for exploratory analysis (maybe this is the case??). These counts as I understand are the residuals of the model, so not really sure what to do.
  
  Another option is to use the normalised counts produced by EdgeR, but these ones doesn't account for the unwanted variation (W\_1, W\_2 covariates) (right???)
  
  Any advice or nice solution?? :=)
  Thanks in advance!!
  
  **Tags:** rna-seq, normalization, edger
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2022-07-18. Score: 2 | Accepted answer
  
  You could correct your logCPMs (`edgeR::cpm(y, log=TRUE)`) with the covariates from RUVseq using limma and then use these for plotting:
  
  ```
  set2 <- RUVg(...)
  
  designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))
  
  y <- DGEList(counts=counts(set2), group=pheno$condition4)
  y <- calcNormFactors(y, method="upperquartile")
  y <- estimateGLMCommonDisp(y, designRUV)
  y <- estimateGLMTagwiseDisp(y, designRUV)
  fit <- glmFit(y, designRUV)
  
  # correction with limma
  logcpm <- cpm(y, log=TRUE)
  design <- model.matrix(~0+condition4, pData(set2))
  corrected <- removeBatchEffect(logcpm, design=design, covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])
  
  ```
  
  > **Comments:**
  > **username_0**: Thank you @username_1. Can I ask why you used a "-" in the `model.matrix`? Also, after reading your answer I have been checking the `removeBatchEffect` documentation and when you put the covariates argument you are setting all the phenoData of set2... Don't would be `covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])` ?. Thank you!!!  
  > **username_1**: was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit.  
  
  ---
  
  