Limma time-course: how to apply differential contrast if time zero is considered both as treatment and control?
Goal: We are using limma to analyze time course experiment with 3 time points, and one treatment vs. control. We would like to test this combination of time course and one cold stress as follows: Amb_0hr (no stress time zero), Amb_1hr, Amb_6hr, Cld_1hr (cold stress after 1 hour), Cld_6hr.
Experiment description: we started with 5 flasks of bacteria; extracted RNA for time zero from one of them; applied cold stress only on two of them; extracted RNA after 1hr from one stressed and one control flasks; extracted RNA after 6 hours form one stressed and one control flasks (everything was done in triplicates).
Issue": In the contrasts we used (below) time zero is considered both as "time zero none stress", and "time zero cold stress". PROBLEM: If we use the differential contrast suggested below [Dif1_hr =( Cld _1hr-Amb_0hr)-(Amb_1hr-Amb_0hr)] we get more upregulated genes then expected in "standard" contrast" (Cld _1hr-Amb_0hr). In specific, we get 204 upregulated genes for the differential contrast Dif1_hr. However using "standard" contrast we get only 79 upregulated genes (see results table below).
A.** This is the differential contrast:**
cont.wt <- makeContrasts(
Dif1_hr =(Cld_1hr-Amb_0hr)-(Amb_1hr-Amb_0hr),
Dif_6hr=(Cld _6hr- Cld _1hr)-(Amb_6hr-Amb_1hr),
levels=design)
B.** This is the "standard" contrast":**
cont.wt <- makeContrasts( " Cld _1hr-Amb_0hr"," Cld _6hr- Cld _1hr", "Amb_1hr-Amb_0hr","Amb_6hr-Amb_1hr",levels=design)
C. Results Table:
Problem: We get 204 upregulated genes for the differential contrast Dif1_hr. However using "standard" contrast (Cld _1hr-Amb_0hr) we get only 79 upregulated genes:
Dif1_hr Dif6_hr Cld_1hr-Amb_0hr Cld_6hr- Cld_1hr Amb_1hr-Amb_0hr Amb_6hr-Amb_1hr
54D/204U 212D/138U 367D/79U 311D/368U 334D/61U 42D/55U
*D=Downregulated, U=upregulated.
What is the best contrats to ask "Which genes respond at either 1 hour or 6 hour of cold stress?
Did I make a correct differential contrast?
Detailed ANALYSIS:
1. Data taken from RNA-seq experiment.
2. Design matrix (two time points :1,6 hr; treatments: control, coldStress, triplicate for each treatment):
Amb_0hr Amb_1hr Amb_6hr Cld_1hr Cld_6hr
1 1 0 0 0 0
2 1 0 0 0 0
3 1 0 0 0 0
4 0 1 0 0 0
5 0 1 0 0 0
6 0 1 0 0 0
7 0 0 1 0 0
8 0 0 1 0 0
9 0 0 1 0 0
10 0 0 0 1 0
11 0 0 0 1 0
12 0 0 0 1 0
13 0 0 0 0 1
14 0 0 0 0 1
15 0 0 0 0 1
3. Commands used:
lev<-c("Amb_0hr","Amb_1hr","Amb_6hr","Cld_1hr","Cld_6hr")
all_groups_repeats<-c("Amb_0hr", "Amb_0hr", "Amb_0hr", "Amb_1hr", "Amb_1hr", "Amb_1hr", "Amb_6hr", "Amb_6hr", "Amb_6hr", "Cld_1hr", " Cld _1hr", " Cld _1hr", " Cld _6hr", " Cld _6hr", " Cld _6hr")
f <- factor(all_groups_repeats, levels=lev)
design <- model.matrix(~0+f)
colnames(design) <- lev
fit <- lmFit(eset, design)
v <- voom(xCurrent, design, plot=FALSE)#Voom fit
fit <- lmFit(v , design)
Contrast used:
cont.wt <- makeContrasts(
Dif1_hr =(Cld _1hr-Amb_0hr)-(Amb_1hr-Amb_0hr),
Dif_6hr=(Cld _6hr-Cld _1hr)-(Amb_6hr-Amb_1hr),
levels=design)
fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
#Final DE analysis
wt<-decideTests(fit2,p.value = 0.05, adjust.method = "BH", method = "global")
summary(wt)#Produce summary table
Thank you very much!!
Unfortunately, our problem is not solved yet. Could you Please suggest how should we find which genes respond to 1hr of cold stress?
Reminder: Experiment description: we started with 5 flasks of bacteria; extracted RNA for time zero from one of them; applied cold stress only on two of them; extracted RNA after 1hr from one stressed and one control flasks; extracted RNA after 6 hours form one stressed and one control flasks (everything was done in triplicates).
Issue": If we want to find which genes respond to 1 hr cold stress we would like to use the "differential contrast" suggested in page 47 of the user guide. While it is straightforward for 6hr cold stress (i.e., makeContrasts( Dif_6hr=(Cld _6hr- Cld _1hr)-(Amb_6hr-Amb_1hr), levels=design) it is not trivial for 1hr cold stress. The reason for that is that, in our time zero is considered both as "time zero none stress", and "time zero cold stress".
Considering this. Coluld you please suggest what is the best way to build differential contrast to 1hr hour cold stress?
I was thinking about the two following options:
1. For 1 hour cold stress:
Use this "standard contrat":
cont.wt <- makeContrasts( "Cld _1hr-Amb_0hr"," Cld _6hr- Cld _1hr", "Amb_1hr-Amb_0hr","Amb_6hr-Amb_1hr",levels=design)
And then to substract the list of over expressed genes found in "Amb_1hr-Amb_0hr" from the list of genes found by "Cld _1hr-Amb_0hr".
Is this a correct approach?
2. For 6 hour cold stress:
I have two options:
a. Using the same appraoch used for 1hr cold stress based on the "standard contrat".
OR
b. Using the differential contrast based on page 47 of the user guide:
cont.wt <- makeContrasts( Dif_6hr=(Cld _6hr- Cld _1hr)-(Amb_6hr-Amb_1hr), levels=design)
Which one is better for 6 hours?