differential gene expression analysis

Question

[ATCG]

I'm trying to do differential gene expression analysis (RNA-Seq data) using Limma and EdgeR How can I adjust my count matrix for technical confounders(RIN, batch)?

a. Is it possible to adjust for confounders within Limma?

b. Can I adjust the count matrix for covariates using linear regression and import the residuals for downstream differential expression analysis using Limma or EdgeR? Thanks!

[ATCG]

"Yes, typically one simply adds batch and RIN as covariates if they seem relevant. " Do you mean that you can add covariates into the model matrix e.g. design<-model.matrix(~group+batch+RIN) Thanks!

[Devon Ryan]

Yes, correct

[ATCG]

Is there a way to visualize the data (residuals) using PCA or MDS plot after performing the linear modeling? I want to do this to confirm the samples no longer cluster by batch.

Answer 1

[Devon Ryan]

1. Yes, typically one simply adds batch and RIN as covariates if they seem relevant. You could also use combat from the SVA package if really necessary. In limma there's also a removeBatchEffects() function, though that wouldn't handle RIN.
2. See the SVA package in R.

[Gordon Smyth]

Actually removeBatchEffects() can handle numeric covariates like RIN easily enough. Just set batch=batch and covariates=RIN in the removeBatchEffects call. The covariates argument wasn't an original part of the removeBatchEffects function, but was introduced with the April 2012 Bioconductor release.

Answer 2

[ATCG]

Devon, How can I plot the data going into the analysis? design<-model.matrix(~group+batch+RIN) will include RIN in the model and this should change the way samples cluster by PCA or MDS plot. Is there a way to extract the corrected counts matrix? Thank you!

[Devon Ryan]

Including things in a model will have absolutely no effect on PCA or MDS. Check some residuals or some scatter plots if you want to QC how well the correction is being done.

[Gordon Smyth]

There's no such thing as a corrected counts matrix, but one can correct the logCPM values. To see visually how the samples cluster after adjustment for batch and RIN:

logCPM <- cpm(y,log=TRUE)
logCPM.corrected <- removeBatchEffect(logCPM, batch=batch, covariates=RIN,
                                  design=model.matrix(~group))
plotMDS(logCPM.corrected)

I would do this the same way regardless of whether you are using limma or edgeR for the DE analysis. Note that logCPM.corrected is only for plotting purposes --- don't use it for the DE analysis.