Hello,
I was using ComBat-seq from the sva package to correct the batch effects for my RNA-seq samples, which were from two different batches but exactly the same platform.
Before doing this, I've checked my data for batch effects, here is the code I used and returned results:
library(DESeq2)
sample_info <- data.frame(sample = colnames(count_mat),
condition = factor(rep(c("pre_treat", "post_treat"),
times = 9)),
batch = factor(c(rep(1, 8), rep(2, 10))))
rownames(sample_info) <- sample_info$sample
dds <- DESeqDataSetFromMatrix(countData = count_mat,
colData = sample_info,
design= ~ condition+batch)
pca_dat <- plotPCA(DESeqTransform(dds), intgroup=c("condition", "batch"),
returnData=TRUE)
I did the batch effects correction with ComBat-seq , and the input counts was the raw counts matrix, similarly, I also plotted the PCA:
gene_exp_adj <- ComBat_seq(counts = counts, batch = sample_info$batch,
group = sample_info$condition)
but apparently, the results didn't change much:
I also performed same analysis on my lncRNA datasets, and the results of before/after batch correction were very similar. Does anyone know why? Or is there another tool I could try? Any advice would be appreciated!
1 answer
First of all you should follow the DESeq2 manual and use plotPCA correctly. You are giving it explicitely a DESeqTransform object (the manual does not suggest that -- it also makes no sense) and the axis limits of the PCA indicate that data are neither log-transformed - and based on the code probably not normalized as well. Please read the manual, do exactly what it does for the PCA and repeat the plots. My guess is that you are doing PCA on raw counts which makes no sense at all. So manual => DESeqDataSet() => vst() => plotPCA().
Thank you for your enthusiastic answer! Well, after normalising the original data, I got the pca results as below:
gene datasets-before batch correction:
gene datasets-after batch correction:
Also, I tried removeBatchEffect function from limma package on gene datasets, but how to determine which one has better batch correction effect?
batch correction for gene datasets with limma:
But combat_seq didn't handle my lncRNA datasets well.
lncRNA datasets-before batch correction
lncRNA datasets-after batch correction
So I also tried removeBatchEffect function from limma package, it did seem to work:
batch correction for lncRNA datasets with limma:
Still not knowing the reason why lncRNA datasets have more obvious batch effect than mRNA datasets, since both of them were sequenced from the exact same source.
And why combat_seq failed to correct batch effects for lncRNA datasets.
Log in to answer this question.
-
Batch Effect Correction in RNA-seq Data Using SVA and DESeq2
written by Nelo •
-
DEseq2 on expected counts
written by lyan125 •
-
How can I conduct a meta-analysis of RNA-seq or microarray data?
written by wyt1995 •
-
DESeq2: When to split multiple group samples for 'more' accurate sizefactors, dispersions etc.
written by kelen
-
DESeq2 Batch correction enquiry?
written by prithvi.mastermind •
-
Correct for batch-effects which are not visible in pca-plot of raw-data (RNA-Seq)
written by Tobias.Wohland
What do you mean when you say it didn't work. Was there an error, or was the batch effect not resolved?
The later case, the batch effects didn't seem to be resolved.
Are you sure there’s any actual batch effects to begin with? I’m not an RNAseq expert but it doesn’t seem like it.
Maybe the results for my lncRNA datasets is more obvious, nothing seems to have changed: before batch correction
after batch correction