Dear all,

I am analysing some data from a DNA methylation array study. The comparison of interest is between DNA methylation in neurons that were differentiated from neuronal precursor cells (NPCs) made by differentiating iPSCs from cases (n = 4) and controls (n = 4). For each individual, between two and six DNA samples were assessed. For a given individual, these DNA samples were obtained from neurons that were grown on separate occasions for varying lengths of time (different passage numbers). For some individuals, methylation was measured in neurons differentiated from independent derivations of NPCs. Unfortunately a lot of the information on the passage number of the neurons and the identity of the NPC source of the neurons was not supplied. 

This means that my dataset looks like this:

I have been trying to work out how best to treat the replicate arrays from each individual. The fact that the replicate arrays from each individual are not strictly technical replicates has made me reluctant to simply average them prior to carrying out analysis of differential methylation (although I have tried this approach). I, therefore, wondered if duplicateCorrelation, with individual as the blocking factor, might be a reasonable thing to do.

So far, I have carried out the case-control comparison in two ways:

1. averaging across the arrays for each individual, as follows (averaged_phenodata is a matrix where each individual is represented by one row and averaged_meth_data is average array data for each individual obtained using avearrays):

design <- model.matrix(~Case_status + Sex, data = averaged_phenodata)

fit <- lmfit(averaged_meth_data, design)

2. Using duplicateCorrelation, as follows:

design <- model.matrix(~Case_status + Sex, data = phenodata)

corfit <- duplicateCorrelation(meth_data, design, block=phenodata$Individual)

fit <- lmFit(meth_data, design, block=phenodata$Individual, correlation=corfit$consensus)

The results obtained from the two approaches are similar in terms of the identities of the top ranked loci but duplicateCorrelation results in many more significant results. 

My questions are:

1. Is this, in theory, a legitimate use of duplicateCorrelation? Does the fact that the expected correlation between the replicates from each individual is not necessarily the same (due to differing levels of similarity of the replicates) invalidate the use of the consensus correlation?

2. Have I gone about implementing it correctly?

Any guidance would be very much appreciated.

Thanks,

Rosie

Yes, the unknown number of passages could be an issue as it may introduce additional correlations between samples. For example, samples from different patients may be correlated if they have been passaged a similar number of times - this would be totally missed by the model. This may or may not be a realistic concern, based on the effect of passaging in your particular system. If the effect is weak, it may be fine to ignore it.

In any case, the use of duplicateCorrelation is best suited to situations where correlations are present across different levels of experimental factors (e.g., sex, group). Here, your are modelling correlations between samples that belong to the same patient, which have the same levels based on the level of the patient.

So, you might as well average the observations for each sample; this avoids the assumptions (and computational work) involved in running duplicateCorrelation. Averaging may also provide a little bit of protection against the passage problem. If you can assume that the number of passages is randomly distributed across samples, it should cancel out between patients if you take the average of samples within each patient.

There are several additional points to note with your first approach. The first is that you can explicitly block on the Slide number, which may be a batch effect contributing to some unwanted variance. The second is that you can run arrayWeights to account for the fact that patients have different numbers of averaged samples (and presumably different precision, i.e., more samples = less variance).

See also A: limma - technical replicates: duplicateCorrelation() or avereps()?.