Question: Time-series Analysis using DESeq2 and LRT test
0
gravatar for Assa Yeroslaviz
2.7 years ago by
Assa Yeroslaviz1.4k
Munich, Germany
Assa Yeroslaviz1.4k wrote:

Hi all 

I have a data set of 27 sample from 9 different time-points with triplicates of control and knockouts. 

sampleNames    condition    treatment    TP    batch
A    CTRL0    CTRL    0    1
B    CTRL0    CTRL    0    1
C    CTRL0    CTRL    0    1
D    KO0    KO    0    1
E    KO0    KO    0    1
F    KO0    KO    0    1
G    CTRL3    CTRL    3    1
H    CTRL3    CTRL    3    1
I    CTRL3    CTRL    3    1
J    KO3    KO    3    1
K    KO3    KO    3    1
L    KO3    KO    3    1
M    CTRL6    CTRL    6    1
N    CTRL6    CTRL    6    1
O    CTRL6    CTRL    6    1
P    KO6    KO    6    1
P2    KO6    KO    6    2
Q    KO6    KO    6    1
Q2    KO6    KO    6    2
R    KO6    KO    6    1
R2    KO6    KO    6    2
S    CTRL9    CTRL    9    1
T    CTRL9    CTRL    9    1
U    CTRL9    CTRL    9    1
V    KO9    KO    9    1
W    KO9    KO    9    1
X    KO9    KO    9    1

I have already done a pair0wise analysis within each time-point, comparing the KO against the control. I would like now to test for genes differentially expressed in over time. 

from previous experiment I know I can use the LRT test for that. This time it is different though, as I have control samples for each of the time points. in the previous analyses I have had a 0 time-point and have therefore used the reduced=~1 reduced model. 

In the case above I would think the correct model will be 

full_model = formula(~ treatment + TP + treatment:TP)

the reduced model I would like to use would be in this case

reduced_model = formula(~ treatment + TP)

And than run :

dds = DESeq(dds, test="LRT", full=full_model, reduced=reduced_model)

Am I correct in my assumption, that this will give me the genes differentially changed in a condition-specific manner over time?

For that I have two more questions - 

1. should I use the condition column rather than the treatment column, where the each control and treated samples are time-point specific?

2. What should be changed, fig i would also like to include possible batch effects? Some samples were sequenced separately, due to technical difficulties.

would it be enough just add the batch effect to the two models like that?

full_model = formula(~ treatment + TP + batch + treatment:TP)

the reduced model I would like to use would be in this case

reduced_model = formula(~ treatment + TP + batch)

thanks

Assa

timecourse deseq2 lrt • 1.9k views
ADD COMMENTlink modified 2.7 years ago • written 2.7 years ago by Assa Yeroslaviz1.4k

I see 4 time points not 9. Can you explain?

You only have 3 samples from batch 2 correct? Do these samples separate in the PCA plot?

ADD REPLYlink written 2.7 years ago by Michael Love25k

1. sorry, yes these are only four time-points. I mixed it with the time-point 9.

2. Yes. the sample from the different batches fit perfectly well together in a PCA

BTW, does it matter what I put in the design parameter, when creating the dds object?. I was thinking of doing it like that:

dds<-DESeqDataSetFromMatrix(countData=countTable, 
                            colData=phenotype, 
                            design= ~ treatment + TP + batch + treatment:TP
)
ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by Assa Yeroslaviz1.4k
Answer: Time-series Analysis using DESeq2 and LRT test
1
gravatar for Michael Love
2.7 years ago by
Michael Love25k
United States
Michael Love25k wrote:

Don't use the 'condition' column above. You just want the 'pure' variables treatment (CTRL and KO) and time (0,3,6,9), not the two pasted together. And make sure time is a factor:

class(dds$time)

Same for batch:

class(dds$batch)

Then you can use the designs you have above with batch, which should help if, as you said the batch 2 samples separate in a PCA.

ADD COMMENTlink written 2.7 years ago by Michael Love25k

thanks for the fast reply.

ADD REPLYlink written 2.7 years ago by Assa Yeroslaviz1.4k
Answer: Time-series Analysis using DESeq2 and LRT test
0
gravatar for Assa Yeroslaviz
2.7 years ago by
Assa Yeroslaviz1.4k
Munich, Germany
Assa Yeroslaviz1.4k wrote:

Hi Michael, 

I am not sure, I have the correct parameters. I have tried as discussed above and run DESeq2 with the following parameters

dds<-DESeqDataSetFromMatrix(countData=countTable, 
                            colData=phenotype, 
                            design= ~ treatment + TP + treatment:TP)
full_model = formula(~ treatment + TP + treatment:TP)
reduced_model = formula(~ treatment + TP)
dds = DESeq(dds, test="LRT", full= full_model, reduced=reduced_model, parallel = TRUE)
res<-results(object = dds, parallel = TRUE, cooksCutoff = FALSE)

But after checking the plotCounts(), it doesn't look like it is really changing over time. I added the two highest significant genes according to the res object.

pic1pic2

 

there is clear time or condition-specific behaviour, unless the subtle changes I can detect are enough for DESeq2 to set these genes as condition-specific changes over time. is this the case?

thanks

assa 

ADD COMMENTlink written 2.7 years ago by Assa Yeroslaviz1.4k

If I look at the data I get the feeling, that the changes are over time, but not condition-specific. Is this the correct approach to identify genes that are significantly differentially expressed between the Knock-out and the control over all time points?

 

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by Assa Yeroslaviz1.4k
The first one looks plausible. It would easier if you colored by WT vs KO. Use the intgroup argument.
ADD REPLYlink written 2.7 years ago by Michael Love25k

What is happening when I am running the results command with any parameters?

these are the resultNames(dds) I have 

> resultsNames(dds)
[1] "Intercept"                "treatment_mad2KO_vs_CTRL"
[3] "TP_3_vs_0"                "TP_6_vs_0"               
[5] "TP_9_vs_0"                "treatmentmad2KO.TP3"     
[7] "treatmentmad2KO.TP6"      "treatmentmad2KO.TP9" 

As far as I understand it, I always get per default the last parameter in the list. So here i get  "treatmentmad2KO.TP9". But is this what i want for the differences between the mad2KO and CTRL over time?Wouldn't this one just give me the differences between TP9 and TP0?

Won't the contrast "treatment_mad2KO_vs_CTRL" be better for what i am looking for?

thanks again

Assa

 

 

ADD REPLYlink written 2.7 years ago by Assa Yeroslaviz1.4k

This is discussed in the man page for ?results, about what results() gives you when you've performed an LRT.

ADD REPLYlink written 2.7 years ago by Michael Love25k
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 16.09
Traffic: 275 users visited in the last hour