# plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates

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!! ## 1 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"])

• 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!!! Jul 19 at 6:46
• was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit. Jul 19 at 12:00