Hi,

I decided to split DE genes over time with DESeq2 into several parts, as I figured out that it is too long. I hope it was a good idea.

I'm doing a time course analysis with deseq2, I have 8 different points and multiple replica.

when I am looking at resultsNames(dds) in my analysis I see

>resultsNames(dds)

# [1] "Intercept"      "replica_2_vs_1" "replica_3_vs_1" "time_16h_vs_0h" "time_1d_vs_0h"
# [6] "time_24h_vs_0h" "time_30h_vs_0h" "time_48h_vs_0h" "time_72h_vs_0h" "time_90h_vs_0h"

but when I do results(dds), I get only the last time point

> res
log2 fold change (MLE): time 90h vs 0h
LRT p-value: '~ replica + time' vs '~ replica'
DataFrame with 17558 rows and 6 columns
...

Does it means that with this coefficient I get all the genes which are changing over all TPs?

How can I interpret the other coefficients? Are the changes over time up to a specific TP (e.g. 16h_vs_0h) means all genes changing from 0h to 16h?

Is it possible to extract with this kind of analysis also all pair-wise comparisons (like 1d vs. 16h, or 24h vs. 48h), or do i need to run a "normal" wald test for the pair-wise comparisons?

the script I run is below:

#reading in the phenotypic data from a file
phenotype <- read.delim2("phenotypeData.txt", colClasses = c(rep("factor", 3)))
# create the DESeq object from a matrix.
dds<-DESeqDataSetFromMatrix(countData=countTable, colData=phenotype, design= ~ replica + time )
dds = DESeq(dds, test="LRT", reduced=~replica)

res<-results(dds)

The results function takes a "contrast" or "name" argument so that you can specify the result you are interested in, cf the help page via ?results

@Steve - thanks for the response. I did read the help page and the user guide, but I'm still not sure if I understand it correctly.

From the list I get from resultsNames() I can use in the results() command with the name= parameter. Everything else I can constract with the contracts= parameter?

Does this means, that the contrast=c("time", "16", "0") and name="time_16h_vs_0h" will give me the same results?

Back from to my question from above. So if I understand it now correctly, I can just use contrast=c("time", "24", "48") to get all genes which are changed between the TP 48 hours to 24 hours, even though it is not in th resultNames() object.

Is it better to use here the Wald test, when comparing only two TP agianst each other?

Assa

Here's the relevant section from ?results

"For analyses using the likelihood ratio test (using nbinomLRT), the p-values are determined solely
by the difference in deviance between the full and reduced model formula. A log2 fold change is
included, which can be controlled using the name argument, or by default this will be the estimated
coefficient for the last element of resultsNames(object)."

There is also a section in the vignette which might help explain things:

"3.5 Likelihood ratio test

DESeq2 offers two kinds of hypothesis tests: the Wald test, where we use the estimated standard error of
a log2 fold change to test if it is equal to zero, and the likelihood ratio test (LRT)..."

For the likelihood ratio test results table, you can just disregard the LFC column and think of this as giving you a single p-value for a set of fold changes.

So for DESeq2:

Wald test for testing specific fold changes or linear combinations (condition specific changes at a particular time, using 'name' and telling results that you want test="Wald")

LRT for testing whether sets of fold changes are all equal to zero (any condition specific changes over time)

You can also take a look at the workflow which has a time series example:

www.bioconductor.org/help/workflows/rnaseqGene/#time

I might suggest you talk to a local statistician if you have further questions on the difference between a Wald test and a likelihood ratio test, or on how to interpret the various results tables. The coefficients given by DESeq2 here are standard for any kind of a linear model.