How to perform batch correction when we have only single condition?
2
2
Entering edit mode
3.8 years ago
svp ▴ 680

I have only single condition "Cancer" and where processed as three batches.

SampleID    Condition   Batch
Sample1      Cancer      A1
Sample2      Cancer      A1
Sample3      Cancer      A1
Sample4      Cancer      A1
Sample5      Cancer      A1
Sample6      Cancer      A1
Sample7      Cancer      A1
Sample8      Cancer      A2
Sample9      Cancer      A2
Sample10     Cancer      A2
Sample11     Cancer      A2
Sample12     Cancer      A2
Sample13     Cancer      A2
Sample14     Cancer      A2
Sample15     Cancer      A3
Sample16     Cancer      A3
Sample17     Cancer      A3
Sample18     Cancer      A3
Sample19     Cancer      A3
Sample20     Cancer      A3
Sample21     Cancer      A3

followed following method to normalize and batch correction

#Read the count table data
cts <- read.table("raw_read_count.txt", header=TRUE, row.names=1)

#Sample information and barch details
coldata <- read.table("sample_info.txt", header=TRUE, row.names=1)

#Create a DESeq2 object named dds from the gene read count and sample information
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ 1)

#assign the batch column
dds$batch <- factor(dds$Batch)

#estimate the library size correction and save the normalized counts matrix
dds <- estimateSizeFactors(dds)
norm.cts <- counts(dds, normalized=TRUE)

#perform regularized log transformation (blind false)
rld_f <- rlog(dds, blind=FALSE)

Next step is to do a batch correction. I followed the method given here as follows:

library(sva)
mm <- model.matrix(~ Condition, colData(dds))
mm0 <- model.matrix(~ 1, colData(dds))
norm.cts <- norm.cts[rowSums(norm.cts) > 0,]
fit <- svaseq(norm.cts, mod=mm, mod0=mm0, n.sv=2)

where it ask for mm and mm0; and mm does not work for me as I have only single condition. How can I go ahead with this? What changes has to made in the above code ?

sva combat batch-effect DESeq2 • 1.0k views
ADD COMMENT
2
Entering edit mode
3.7 years ago
svp ▴ 680

I am updating the answer. If anyone stuck on the operation can follow

#get the matrix of regularized log count
rlog_count = assay(rld_f)

#Create model
cb.corr.model <- model.matrix(~1, data = coldata)

#Batch correction using combat
cb.corr.counts = ComBat(dat=rlog_count,
                        batch=dds$batch,
                        mod=cb.corr.model,
                        par.prior=TRUE,
                        prior.plot=FALSE)
ADD COMMENT
1
Entering edit mode
3.8 years ago

SVA is used to find (not necessarily remove) unknown sources of variation. Here your source of error is known - you know which batch each sample came from.

I recommend using removeBatchEffects from limma instead.

ADD COMMENT

Login before adding your answer.

Traffic: 2981 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6