I have the following experimental design and want to estimate the differences between nonallergic.BPE vs nonallergic.Medium, allergic.BPE vs. allergic.Medium etc. I would like to use patient as a blocking factor
patient genotype condition 73644 1 nonallergic Medium 73645 1 nonallergic BPE250 73646 2 nonallergic Medium 73647 2 nonallergic BPE250 73648 3 nonallergic Medium 73649 3 nonallergic BPE250 73650 4 allergic Medium 73651 4 allergic BPE250 73652 5 allergic Medium 73653 5 allergic BPE250 73654 6 allergic Medium 73655 6 allergic BPE250 73656 7 allergic Medium 73657 7 allergic BPE250 73877 8 nonallergic Medium 73878 8 nonallergic BPE250
In limma-voom it did it like this:
group <- factor(paste(genotype, condition, sep=".")) design <- model.matrix(~0 + group) vm <- voomWithQualityWeights(y, design=design) corfit <- duplicateCorrelation(vm, design, block = ExpDesign$patient) vm <- voomWithQualityWeights(y, design, block = ExpDesign$patient, correlation = corfit$consensus, plot= TRUE) vfit <- lmFit(vm, design, block=ExpDesign$patient,correlation=corfit$consensus) my.contrasts <- makeContrasts( aMvsnaM = allergic.Medium-nonallergic.Medium, aBPEvsnaBPE = allergic.BPE250-nonallergic.BPE250, aBPEvsaM = allergic.BPE250-allergic.Medium, naBPEvsnaM = nonallergic.BPE250-nonallergic.Medium, levels=design)
When I read about the design formula and matrix in edgeR, I understood that the blocking should be set up like this:
design <- model.matrix(~0+ patient + group)
Unfortunately this gives me a "Design matrix not of full rank" error for following matrix:
patient1 patient2 patient3 patient4 patient5 patient6 patient7 patient8 groupallergic.Medium groupnonallergic.BPE250 groupnonallergic.Medium 1 1 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 0 3 0 1 0 0 0 0 0 0 0 0 1 4 0 1 0 0 0 0 0 0 0 1 0 5 0 0 1 0 0 0 0 0 0 0 1 6 0 0 1 0 0 0 0 0 0 1 0 7 0 0 0 1 0 0 0 0 1 0 0 8 0 0 0 1 0 0 0 0 0 0 0 9 0 0 0 0 1 0 0 0 1 0 0 10 0 0 0 0 1 0 0 0 0 0 0 11 0 0 0 0 0 1 0 0 1 0 0 12 0 0 0 0 0 1 0 0 0 0 0 13 0 0 0 0 0 0 1 0 1 0 0 14 0 0 0 0 0 0 1 0 0 0 0 15 0 0 0 0 0 0 0 1 0 0 1 16 0 0 0 0 0 0 0 1 0 1 0 attr(,"assign") [1] 1 1 1 1 1 1 1 1 2 2 2 attr(,"contrasts") attr(,"contrasts")$patient [1] "contr.treatment" attr(,"contrasts")$group [1] "contr.treatment"
How do I have to set up my matrix in order to get the differences for my contrasts?
Unfortunately I also don't know whether I should use a model with intercept or not, as we don't have any biological data that would suggest/dissuade that the treatment would have more effect in allergic vs nonallergic patients... Do you have any advice on it? How could I figure out which model would be the right one?
Cheers
Thank you for your explanation Aaron and the appropriate section in the guide. I don't have to use edgeR, but as we don't know what outcome we expect, I thought it might be a good idea to do the analysis with three packages (edgeR, DeSeq2 and limma-voom) to compare the resulting genes.
I should be able to compare the subseted data from edgeR with the analysis from limma-voom, right?
Is the analysis/code in limma-voom correct? I was just wondering because I don't get any upregulated interleukin genes (which I somewhat expected)...
With this setup don't get a common dispersion factor. is that to be expected?
I'll answer your questions, but first, a rant.
I'm often bemused by people who want to compare different DE analysis methods on data sets where they do not know the ground truth. If the methods yield different DE lists, it's impossible to say which method is more correct. In fact, it's possible for two DE lists to both be correct (i.e., controls the FDR below the nominal threshold) yet also very different in size and composition, depending on each method's power for particular types of genes, e.g., low/high-abundance, low/high-variance. One hopes that the differences would decrease as the sample size increases, but most RNA-seq experiments are underpowered and different methods will use different assumptions to "fill in the gaps", so to speak.
The justification I often hear for using multiple methods is something like "if we take the intersection of DE lists from all methods, we'll get a confident set of DE genes". This is a specious argument. If I want a more confident set of DE genes, I lower the FDR threshold and trust that my method is controlling the FDR correctly. If it isn't, then no amount of intersections will help. In fact, a naive intersection between DE lists (each detected at a certain FDR threshold) is not guaranteed to yield a set of genes where the FDR is controlled below that same threshold! To be correct, we need to use intersection-union tests instead, and I'm pretty sure that no one who advocates for intersections is actually doing that. Finally, an intersection of any kind will reduce power; a gene that fails to be detected in any method (e.g., because a particular method is underpowered for a particular expression profile) will not be in the intersection, even if it is a true positive.
Okay, that's enough ranting. Onto your questions.
I should be able to compare the subseted data from edgeR with the analysis from limma-voom, right?
Yes and no. Yes, because it is the same contrast. No, because the analysis on the subset will naturally have less power than the analysis of the full data set where more information is available for variance/dispersion estimation. So you'll naturally get different results, without even considering the differences in the models.
Is the analysis/code in limma-voom correct?
I would do a second
duplicateCorrelationafter the secondvoom, but otherwise it looks fine to me (assuming you rancontrasts.fitandeBayescorrectly). If expected DE genes are not being detected, I would suggest looking at their expression values directly to see whether it is a modelling issue or if there wasn't a difference in the first place. Your design is still simple enough for direct examination of the log-CPMs.With this setup don't get a common dispersion factor. is that to be expected?
That's because you saved the output of
estimateDisptoy_edge, noty.I'll answer your questions, but first, a rant.
Fair enough and you have a point. I just figured that if the same genes were detected by multiple methods, they were more likely to be identified correctly.
I would suggest looking at their expression values directly to see whether it is a modelling issue or if there wasn't a difference in the first place
I've had a look at the
lcpmfrom a few interleuking genes and the values between the 16 samples differ from e.g 7 to 5. Is that difference not large enough to be picked up?Depends on the variance within each group, compared to the log-fold change between groups.
So I ran the modeling with a second
duplicateCorrelationas suggested, and got acorfit$consensusof0.6093099. I was happy with the results, till I checked the pvalue distributionand the
vm$EDoes that mean the
duplicateCorrelationdidn't work, or do I have more batch effects in my data I haven't removed yet? The MDS plots before and after voom look almost identicalThe MDS plots won't directly respond to any blocking within
duplicateCorrelation. In fact, it's impossible for them to do so, as you can't explicitly estimate and remove the batch effects here without also removing the effect of interest. Any difference will be subtle and due to the effect ofduplicateCorrelationblocking on the precision of the fitted linear model when correlations between samples are taken into consideration.I don't understand what you intend to show with the p-value distribution. If you have DE genes in that contrast, obviously you'll have more p-values below the threshold than expected under the null.