Question: Limma time-course: how to apply differential contrast if time zero is considered both as treatment and control?
3
gravatar for harelarik
17 months ago by
harelarik30
harelarik30 wrote:

Summary: To find which genes respond to 1 hr cold stress we would like to use the "differential contrast" suggested in page 47 of the limma user guide. However it is not trivial, as in our experiment time zero is considered both as "time zero none stress", and "time zero cold stress".**

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": 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 experiment 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?


Potential solutions:

A. For 1 hour cold stress:

I was thinking of using this contrast:

cont.wt <- makeContrasts( "Cld _1hr-Amb_0hr", "Amb_1hr-Amb_0hr", "Cld _6hr- Cld _1hr", "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" (same is done for under expressed).

Is this a correct approach?

B. For 6 hour cold stress:

I have two options:

a. Using the same contrast used for 1hr cold stress above, and then to subsutract the list of over expressed genes found by "Amb_6hr-Amb_1hr" from the list of overexpressed genes found by "Cld _6hr- Cld _1hr".

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?



PLEASE AVDICE:

WHAT IS THE BEST CONTRATS TO ASK "WHICH GENES RESPOND AT EITHER 1 HOUR OR 6 HOUR OF COLD STRESS?

Thank you.



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

  1. 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

limma time course • 908 views
ADD COMMENTlink modified 16 months ago by Kevin Blighe41k • written 17 months ago by harelarik30
0
gravatar for Kevin Blighe
16 months ago by
Kevin Blighe41k
Guy's Hospital, London
Kevin Blighe41k wrote:

Hi harelarik,

From my experience of using limma, it's best to keep the model formulae as simple as possible. Well, nothing against limma, really, as an overly complex formula will produce unreliable information in any regression analysis. One of my main worries in your situation is the very low sample numbers.

Given the relatively complex set-up and the low sample numbers, I believe that you will find it extremely difficult to arrive at just a single list of statistics and be able to say: 'This is my solution for 6 hour cold stress', et cetera.

I would begin by critically thinking about what each comparison is actually showing. For example, Cld _6hr-Cld _1hr will show me genes that are statistically differentially expressed between the 1 hour and 6 hour time-points in the cold-stress samples, whereas Amb_6hr-Amb_1hr will show me the same but at ambient / room temperature. If I then compare these two listings (even by eye), I will gain insights into genes that are activated or repressed when cold-stress is applied.

Then do the same for 1 and 0 hours, and 6 and 0 hours. You will end up with 6 different gene listings, each showing various things.

--------------------------------

By the way, in addition and considering that you have bacterial RNA-seq, you may be interested in taking a look at Rockhopper. It runs quickly and will return FASTA sequences that are statistically differentially expressed (and provide P and Q values). You can then infer functionality of these sequences via blastx.

Kevin

ADD COMMENTlink modified 16 months ago • written 16 months ago by Kevin Blighe41k

Thank you very much, Kevin! We have 3 biological repeats per samples so: 21 samples.

Cld _6hr-Cld _1hr will show me genes that are statistically differentially expressed between the 1 hour and 6 hour time-points in the >cold-stress samples, whereas Amb_6hr-Amb_1hr will show me the same but at ambient / room temperature. If I then compare >these two listings (even by eye), I will gain insights into genes that are activated or repressed when cold-stress is applied.

That was my initial thought. But then I saw the "differential contrast" suggested in page 47 of the limma user guide, and I thought it might have some advantage. I.e., is using: Dif_6hr=(Cld _6hr- Cld _1hr)-(Amb_6hr-Amb_1hr), better then what you suggested above?? The approach you have suggeted, sounds like the same thing, but I do not get the same number of differentially experssed genes in the final result.

ADD REPLYlink modified 16 months ago • written 16 months ago by harelarik30

Yes, you will not get the exact same results because the model comparisons are different. The way that I suggested just produces statistics for each pairwise comparison, whereas your model, Dif_6hr=(Cld _6hr- Cld _1hr)-(Amb_6hr-Amb_1hr), is attempting to do a further comparison within the model. I just feel that this is too complex for the low numbers of samples that you have, and will produce unreliable results.

The way that I and others with whom I've worked have done these types of studies is through just the simple pairwise comparisons, looking over the results, and then piecing together a storyline based on these.

ADD REPLYlink written 16 months ago by Kevin Blighe41k

Very interesting question @harelarik, I'm facing the same problem right now. In case you had enough samples, would it be correct the makeContrasts( Dif_1hr=(Cld _1hr- Amb_0hr)-(Amb_1hr-Amb_0hr), levels=design) approach for checking the genes implicated in the 1 hour cold stress response?

ADD REPLYlink modified 9 months ago • written 12 months ago by mfalco0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1390 users visited in the last hour