Question: How can Differentially Expressed Genes be found by limma package after removing batch effect by ComBat function in sva package?
0
gravatar for amirmehrgou
12 days ago by
amirmehrgou0 wrote:

I wish to raise a question regarding the detection of differentially expressed genes from combined datasets. So as to merge two expression microarray datasets, following the combination of series matrices and their corresponding annotation files, I used the ComBat function in sva package to remove the batch effect. However, as a beginner bioinformatician, I have almost no idea how to take advantage of the retrieved output of ComBat function of sva package in the Limma package in order to find differentially expressed genes. Therefore, some errors occur usually while I am trying to do so. Accordingly, I pasted an excerpt from my codes in the following to gain your precious help in providing me easy-to-do R codes for approaching my purpose.

If you have time to guide me I will be immensely grateful.

I am looking forward to your response.

Sincerely yours,

Amir Mehrgou

*An excerpt from my code:

library(sva)
library(limma)

allc <- ComBat(as.matrix(all), batch = batch)
allm <- allc - rowMeans(allc)
allm$description <- factor(gr)
design <- model.matrix(~ description + 0, as.data.frame(allm))
colnames(design) <- levels(factor(gr))
fit <- lmFit(as.matrix(allm), design)
cont.matrix <- makeContrasts(Case-Control, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2, 0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="B", number=Inf)
batch effect limma sva R • 99 views
ADD COMMENTlink modified 12 days ago by RamRS27k • written 12 days ago by amirmehrgou0
1

Regardless of the tool used to remove batch effect, I believe the general recommendation is to incorporate the batch as a covariate in the model, and only batch-adjust your matrix for data plotting/exploration.

ADD REPLYlink written 11 days ago by Papyrus330

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.
code_formatting

ADD REPLYlink written 12 days ago by RamRS27k

I was not aware of this method of code writing. Thank you for your guidance.

ADD REPLYlink written 11 days ago by amirmehrgou0

Why are you using ComBat when limma has a removeBatchEffect function?

ADD REPLYlink written 12 days ago by RamRS27k

As you suggested I utilized removeBatchEffect function in limma package. However, at the time of finding differentially expressed genes I encountered the same errors as I had observed when sva package was used. My codes and errors are as follow:

> allB <- removeBatchEffect(as.matrix(all), batch)

> allB$description <- factor(gr)

Warning message:
In allB$description <- factor(gr) : Coercing LHS to a list

> design <- model.matrix(~ description + 0, as.data.frame(allB))

> colnames(design) <- levels(factor(gr))

> fit <- lmFit(as.matrix(allB), design)

Error in rowMeans(y$exprs, na.rm = TRUE) : 'x' must be numeric

> cont.matrix <- makeContrasts(Case-Control, levels=design)

> fit2 <- contrasts.fit(fit, cont.matrix)

Error in contrasts.fit(fit, cont.matrix) : 
  Number of rows of contrast matrix must match number of coefficients in fit
In addition: Warning message:
In contrasts.fit(fit, cont.matrix) :
  row names of contrasts don't match col names of coefficients

> fit2 <- eBayes(fit2, 0.01)

> tT <- topTable(fit2, adjust="fdr", sort.by="B", number=Inf)

Additionally, I made a comparison between removing batch effect by limma (Left side) and sva (Right side) packages the plot of which is attached below. Comparison between limma and sva package in removing batch effect

Which one is more powerful do you think?

Sincerely yours,

Amir

ADD REPLYlink modified 11 days ago • written 11 days ago by amirmehrgou0

Like Papyrus says, account for batch as a covariate. Removing batch effects is not a good idea, especially if the experiment design is unbalanced.

ADD REPLYlink written 11 days ago by RamRS27k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1792 users visited in the last hour