Search
Question: limma- time course multiple treatments
0
11 months ago by
harelarik0 wrote:

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?

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

The series of commands below correspond to chapter 9.6.1 of limma guide.

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?

modified 9 months ago • written 11 months ago by harelarik0
0
9 months ago by
Aaron Lun21k
Cambridge, United Kingdom
Aaron Lun21k wrote:

Your last cont_1hr in makeContrasts should be cont_6hr, presumably.

Section 9.6.1 of the limma user's guide doesn't use treat or method="nested", so I don't know where you're getting those from. This is fine, but remember that treat will only do pairwise comparisons. If you want to test for any effect of treatment at one or all time points, you will need to use eBayes.

P.S. Your tag doesn't make any sense, this is a limma question.

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 treat instead of eBayes, and to leave lfc at the default value when using decideTests."

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 use eBayes to do so. But it seems like you don't want that, so that's okay.