Dear Bioconductor community,

I encountered a difficult problem in correcting for batch effect in repeated measurement data, which has troubled me for a long time. This is not an experimental design usually discussed on the forum, so I'm having trouble finding solutions.

I am using cohort data, which follows the same group of healthy people and collects plasma samples every 3 years. We have performed small RNA sequencing on 1500+ plasma samples collected from multiple visits, with many of them repeated on the same individuals.

For this particular study, I want to address the following research aims:

1. Assess the long-term variability of miRNA levels in human plasma
2. Investigate whether the variability can be explained by different factors, as well as whether it can predict future development of health conditions (e.g., CHD) using the cohort design

As you can see, this is not a canonical group-level comparison study, but rather I intend to evaluate the within-individual between-visit variability of different miRNAs. I am thinking about using change rate (i.e., the difference in miR divided by years elapsed between 2 visits) to assess this 'variability'.

However, there seems to be some batch effect as shown in the PCA plot below:

Here comes the problem: I found that samples from different visits are not evenly distributed across sequencing batches (figure shown below), so if I use batch-corrected miR data (e.g., with ComBat), I might have removed the variability between visits; if not, I might have batch effect confounded with visits.

I have also tried limma, but it seems that limma does not provide a solution to my problem, which is: how can I correct for batch effect while preserving the variability between visits (or between repeated measures)?

Thank you for your time and patience in reading this post. Any suggestions would be greatly appreciated and invaluable to me.

Best regards,

Wending

There is no free ride in statistics, and since your batch and subject are (mostly) confounded, you will have to make assumptions that may or may not be warranted. If there are any spike-ins used during the processing, you could hypothetically use RUVg from the RUVSeq package to estimate one or more covariates you could add to your design matrix to adjust for batch effects (for example, my core often uses Qiagen's small RNA analysis platform, and there are a bunch of ExiSEQSpikeIn controls that you might be able to use).

But you will never know if you have preserved the variability that you care about while correcting for the batch effect. Ideally the samples would have been randomized appropriately to eliminate the problem, but it's too late for that now.