Batch effect correction using a subset of samples - using DE genes

Question

I have some RNA-seq data with two very obvious batches as you can see in the PCA:

The samples of interest (A - H) are from tumor tissues. In addition, there are data for two cell lines I and J. The cell lines should be the same across the batches, while for the tumor samples there may be some additional biological difference.

My question is, if there is a way to use the cell line data to estimate the technical part of the batch information since there should very little biological difference here. Then I would like to use this to correct the technical aspect of the batch effect in the tumor samples.

When I run normal batch correction like ComBat or add the batch to the DE model the batch will be estimated over all samples including biological differences.

Update

I figured I could determine which genes are differentially expressed between the batches in the cell lines only. As mentioned before, the cell lines are highly similar biologically. Thus these DE genes should represent the technical aspect of the difference.

The plot below shows the PCA when after removing genes DE in the cell lines. Obviously, the differences for the cell lines are greatly reduced. Also the tumor samples are now less clearly separated.

Most of the removed genes were from shorter RNA species and a look at sequencing QC later confirmed that in fact the fragment size in the second batch was lower. Therefore I'm positive that I could identify and remove the technical component of the batch effect.

tl;dr

My data contains a strong difference between batches even for samples that are biologically the same (untreated cell lines). I looked at differentially expressed genes in the cell lines. This revealed to me that the second batch contained a lot more short RNAs due to lower fragment size of the library.

In the end I'm still wondering if there is an elegant solution for this requiring less footwork

Answer 1

The batch effect is expected to be the same across all samples - therefore you want to estimate its effect across all of the samples. Just add batch as an additional parameter to the linear model (this is presuming samples sequenced in the two batches have been randomised).

You can show the effect of "batch effect" removal with a PCA plot by plotting the residual error, e.g. in R if using DESeq:

mypca <-  prcomp(t(assay(varianceStabilizingTransformation(dds)))
pc.resid <- resid(aov(mypca$x~batch.info))
ggplot(pc.resid,aes(x=PC1,y=PC2) + geom_point() 

(I haven't checked the code, but it will be something like this).