Hi,

I'm currently exploring methods for batch correction with my RNA-seq dataset from a series of clinical samples. The goal of this approach is to remove batch effects for PCA and other visualization approaches, as well as to generate input for other analytical pipelines that do not include internal methods for batch correction.

Our data is comprised of time-courses from individual patients that we did not want to separate across batches in order to avoid batch effects when examining within-patient gene expression over time.  We knew this design would make comparing samples across patients challenging, so we included 3 samples that we sequenced as part of every batch (each from a single time point of a different patient).

When examining just the 3 samples repeated across our 5 batches, we see patients separating predominantly by PC1 and batch separating on PC2 . The "removeBatchEffect()" command in limma appears to correct for batch effect very well, with the 5 samples from each patient now clustering very close together and no separation by batch detected by the first 10 PCs.

What I would like to do is use limma batch correction to correct all samples but using only the 3 repeated samples to model the batch correction (since there are certain to be patient to patient differences between batches).

1. Based on Aaron Lun's answer to a similar question (https://support.bioconductor.org/p/91704/), I believe that samples that are contained within only a single batch (all but my repeats) will be used only to estimate variance within that batch, but not used for batch correction. Therefore, RemoveBatchEffects() should only consider these samples when batch correcting.  Is this correct?
2. If my above assumption is incorrect, how would I execute this approach within limma?  Or using another batch correction algorithm (combat, sva, other)?

Thanks for any input the Bioconductor community can provide,

Mike

This is a tricky experimental design. You're basically stuck between two extremes:

• The technical replicates are perfect and have negligible variance, in which case you can manually correct across batches and discard the controls before running limma on the corrected values.
• Technical variance dominates and the contribution of biological variation is negligible, in which case you might as well treat the three controls in each batch as samples in their own right when using limma.

Unfortunately, the truth is probably somewhere in between those two poles, and I'm not aware of any functionality in limma that can smoothly transition from one to the other for a general experimental design. Well, there is, but it's used up in handling the other levels of replication (across patients, and across batches) in your experiment.

The most obvious way to make limma work with off-the-shelf methods is to discard some replication. In particular, you can subset your data so that each batch consists of samples from only one patient - for simplicity, I'll further assume that all samples from one patient belong to one batch. You can then calculate corrected values and discard the controls; fit a linear model to the corrected values using a design matrix containing your factors of interest (e.g., treatment or something, I guess - you didn't specify it); and finally block on the patient effect with duplicateCorrelation. This works because any uncertainty from the batch correction is absorbed into the patient effect and handled by the blocking, given that the patient and batch factors are now equivalent in this subset.

Of course, this is rather unsatisfying if you have to discard a lot of data. I am less sure of other solutions, which would involve weight calculations to account for differences in variance between the controls and other samples. The weights can be interpreted as the reciprocal of the variances of the observations - thus, one could imagine assigning higher weight to control samples that are less variable. Determining the exact value of the weight is the tricky bit. Perhaps arrayWeights might work, but I would have to think it through.

I'm not sure my answer will be helpful since I don't quite follow the variance and weight issues raised by Aaron and may be missing some important aspects here. Take it for what it's worth...

You can try a developmental version of WGCNA from here (dropbox link). To install it, you will need to be able to compile a package with C/C++ code. The package contains function empiricalBayesLM which can also remove batch effects (using approach similar to ComBat). Importantly for your application, the function allows you to specify reference samples (via argument fitToSamples) on which the correction (i.e., the model coefficients) are based. I'd run the function on all your data but specify the "fit to" samples to be the ones that are repeated. You can use either plain linear model adjustment or the empirical Bayes-moderated one which tends to be a bit more stable but leaves a residual batch dependence in the data. Since you sequenced the same 3 samples in each batch, you should not need any "preserved" covariates.

Unlike ComBat, empiricalBayesLM only adjusts means; it does not adjust variances

You should run the function on log-transformed normalized data. The output should be usable for downstream analyses, although it won't help you much with DE analysis in limma.