text: user8886: 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!!
  
  user15793: 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"])
  
  
original_text: 8886: <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>
   15793: <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>
  