Hi,
We would like to test: control, treatment1, and treatment2 in a time course experiment using only few time points.
However, chapter 9.6.1 of limma guide that discusses time-course for few time points, illustrates how to test only two treatments (WT vs. a mutant).
Questions:
1. I wonder if the described approach in this chapter is true also for multiple treatments?
2. If the answer is yes please see what I wrote in section 1-2 below. Is that the correct way to compare multiple treatments over a time course experiment?
Is the contrast.matrix I have used below correct?
-----------------------------------------------------------------
DETAILS:
0. Data from RNA-seq experiment.
1. Design matrix (two time points :1,6 hr; treatments: control, coldStress, heatStress; triplicate for each treatment):
cont_1hr cont_6hr cold_1hr cold_6hr heat_1hr heat_6hr
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 0 0
4 0 1 0 0 0 0
5 0 1 0 0 0 0
6 0 1 0 0 0 0
7 0 0 1 0 0 0
8 0 0 1 0 0 0
9 0 0 1 0 0 0
10 0 0 0 1 0 0
11 0 0 0 1 0 0
12 0 0 0 1 0 0
13 0 0 0 0 1 0
14 0 0 0 0 1 0
15 0 0 0 0 1 0
16 0 0 0 0 0 1
17 0 0 0 0 0 1
18 0 0 0 0 0 1
2. Please see below, Is that a correct way to compare multiple treatments in a time course?
The series of commands below correspond to chapter 9.6.1 of limma guide.
Please note the contrast.matrix.
lev<-c("cont_1hr","cont_6hr","heat_1hr","heat_6hr","cold_1hr","cold_6hr")
f <- factor(all_groups_repeats, levels=lev)
design <- model.matrix(~0+f)
colnames(design) <- lev
v <- voom(rawCounts, design, plot=FALSE)
fit <- lmFit(v , design)
#Is this a correct contrast matrix?
contrast.matrix <- makeContrasts(heat_1hr-cont_1hr,heat_6hr- cont_6hr, cold_1hr- cont_1hr, cold_6hr- cont_1hr, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
fit2<-treat(fit2,lfc=1)#Check parameters!!!
wt<-decideTests(fit2,p.value = 0.1, adjust.method = "BH", method = "nested")
summary(wt)#Produce summary table
Thank you very much for your support!!!
Help anyone?
Thanks a lot!!
* I understand that you approve that it is correct to test multiple treatments using this approach.
You are right - should have been cold_6hr- cont_6hr , Just a mistake.
I meant nestedF (again a mistake). However, it is interesting that limma accepts both nested and nestedF and gives exactly the same results.
I changed the Tags.
Thanks again.
R functions allow partial matching of string arguments, so
"nested"also works. Though fully specifying"nestedF"is safer, because it makes it clearer what you actually want.Thanks!
Regarding "Treat".
Thank you for this comment as well.
I would like to use decideTests() [as it can use nestedF] for final inspection of the results, and to restrict lfc.
One can can use either treat or eBayes, If I am using the latter I need to incorporate lfc in decideests() which could be a problem. As suggested in the decideTests() help: "Although this function enables users to set p-value and lfc cutoffs simultaneously, this is not generally recommended. If the fold changes and p-values are not highly correlated, then the use of a fold change cutoff can increase the false discovery rate above the nominal level. Users wanting to use fold change thresholding are recommended to use
treatinstead ofeBayes, and to leavelfcat the default value when usingdecideTests."This comment is also true for topTable()
That's fine, there is no problem with using
treat. I am just saying that if you wanted to do an ANOVA for any differences between treatments, you will have to useeBayesto do so. But it seems like you don't want that, so that's okay.THANK YOU!