16 months ago by
Cambridge, United Kingdom
I will start by ignoring your questions, because I think you will find it much easier to use:
b <- factor(block)
design <- model.matrix(~b + in.phase + out.phase)
This is permissible as you're expecting oscillations in expression over time, which means that
block is not completely confounded with the time points. A unique solution to the linear model can thus be found using information within each culture. You can then test for significant oscillations by dropping the last two coefficients.
By comparison, the accuracy of
duplicateCorrelation's correlation estimate is dependent on the number of levels of the blocking factor. To illustrate, consider an experiment with 2 batches of 10 samples in each condition:
treatment <- factor(rep(1:2, each=20))
batch <- factor(rep(1:4, each=10))
design <- model.matrix(~treatment)
# Data generation with a batch effect.
dummy <- matrix(rnorm(10000 * length(treatment)), ncol=length(batch))
batch.effect <- matrix(rnorm(10000 * 4), ncol=4)
dummy <- dummy + batch.effect[,as.integer(batch)]
# True value should be 0.5, as both 'dummy' and 'batch.effect' where N(0, 1).
duplicateCorrelation(dummy, design, block=batch)$consensus
You only have two levels in your blocking factor, so you can expect the estimation of the correlation to be similarly inaccurate. This will manifest as anticonservativeness as correlations between samples are not fully considered.
So, long story short, make life easier for yourself and put
block in your design matrix instead.
Now, onto your questions, most of which no longer matter.
3) In practice, yes, due to imprecise estimation of a small positive (or zero) correlation. In theory, probably not; I can't easily imagine how batch effects would manifest as true negative correlations, and it's not surprising that it yields errors downstream. If you get negative estimates of the correlation, it's either because you don't have enough information to estimate the correlation precisely; or the true correlation is very close to zero, such that you don't need to use
duplicateCorrelation at all. It's hard to distinguish between these two possibilities when you only have two levels in your blocking factor. So, put
block in the design matrix instead.
4) Your strategy is flawed for various reasons, but mostly because the null hypothesis is not true for most label permutations. I'll illustrate with a simple two-group example. Imagine I have two groups (A and B) of 10 samples, and I want to shuffle them into mixed groups to test whether my hypothesis tests control the type I error rate. This is only valid for permutations where there are equal numbers of samples from each original group in each mixed group. In any other permutation, the null hypothesis is not guaranteed to be true. For example, if I shuffled such that one group had 9 A's and 1 B, of course I would detect significant differences from the other group containing 9 B's and 1 A (assuming that there were DE genes between A and B in the first place). This means that significant genes cannot be interpreted as false positives.
Another reason is that the variances are inflated when you mix samples from different groups together and pretend that they are "replicates". In and of itself, this is not too much of an issue as the variances are only inflated for genes that are DE between groups - not great, but it should be mostly conservative and should not do too much harm. However, empirical Bayes shrinkage means that information is shared across different genes, which means that the inflated variances can increase the variance of other genes (and conversely, are themselves decreased). This can have unpredictable effects and give a misleading idea of the performance of the tests on real data.
I would suggest you read https://academic.oup.com/nar/article/43/7/e47/2414268, which discusses the limitations of permutation-based methods for hypothesis testing in 'omics contexts.