<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>

CC BY-SA 4.0
2022-07-19T06:46:05.523Z
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!!!
FrAoJm
CC BY-SA 4.0
2022-07-19T12:00:08.643Z
was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit.
alderaan_shot_first
CC BY-SA 4.0
2022-07-18T22:34:44.660Z
2022-07-19T18:19:41.650Z
2022-07-19T18:19:41.650Z
alderaan_shot_first
<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>

<rna-seq><normalization><edger>
plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
# 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 @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!!!  
> **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.  

---

