Extracting DEseq2 results ?resultnames

Question

[Rahul]

Dear Community, I am using the DEseq2 package to analyze my time series data, and referring this vignet (http://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#time-course-experiments)

In brief, the dataset has two groups OP50 (control group) and SC20 (test group), harvested at 7 different time points. we are interested in DE genes in SC20 compared to OP50 at a particular day. for example when comparing SC20_DAY05, I would like to use OP50_Day05 as the reference group. How can I modify the design in DEseq2 to achieve this?

Here Is more information:

1. The metadata

2. My design:

Instead of using interaction term (design= ~Bacteria+Day+Bacteria:Day)

I clubbed variable in the group column and used this:

dds <- DESeqDataSetFromMatrix(countData = countdata, colData = metadata, 
                              design = ~ group)

ddsTC1 <- DESeq(dds, test="LRT", reduced = ~ 1)

3. Here is the output from resultsNames

[1] "Intercept"                      "group_OP50_DAY01_vs_OP50_Day00"
 [3] "group_OP50_DAY03_vs_OP50_Day00" "group_OP50_DAY05_vs_OP50_Day00"
 [5] "group_OP50_DAY07_vs_OP50_Day00" "group_OP50_DAY09_vs_OP50_Day00"
 [7] "group_OP50_DAY11_vs_OP50_Day00" "group_SC20_Day00_vs_OP50_Day00"
 [9] "group_SC20_DAY01_vs_OP50_Day00" "group_SC20_DAY03_vs_OP50_Day00"
[11] "group_SC20_DAY05_vs_OP50_Day00" "group_SC20_DAY07_vs_OP50_Day00"
[13] "group_SC20_DAY09_vs_OP50_Day00" "group_SC20_DAY11_vs_OP50_Day00"

4. Expected outcome

[1] "Intercept"                                      "group_OP50_DAY01_vs_SC20_DAY01"
 [3] "group_OP50_DAY03_vs_SC20_DAY03" "group_OP50_DAY05_vs_SC20_DAY05"
 [5] "group_OP50_DAY07_vs_SC20_DAY07" "group_OP50_DAY09_vs_SC20_DAY09".........

Thanks in advance!

Answer 1

[MB]

I would say the reference condition (OP50_DAY05) should be at the third position:

res <- results(ddsTC1, contrast=c("group", "SC20_DAY05", "OP50_DAY05"))

Alternatively, you should be able to call the coefficient directly via the name argument and geht the same results:

res <- results(ddsTC1, name = "group_SC20_DAY05_vs_OP50_Day00")

Further, I think you should also include the alpha threshold in you hypothesis testing via alpha=. You should check if this applies to LRT, I don't know it from scratch, but normally your default value would be 0.1. If you want to adjust that, adjust that also in the results function.

All the best.

[Rahul]

Thanks for the suggestions. I guess, here,

res <- results(ddsTC1, contrast=c("group", "SC20_DAY05", "OP50_DAY05"))

DESeq2 considers SC20_DAY05 as a reference.

[MB]

No. third position in the vector is the reference. you can also check it by ?results

Answer 2

[Rahul]

I guess I figured it out. But anyone please correct me, If I did something wrong:

I used contrast function

"Intercept" "group_OP50_DAY01_vs_OP50_Day00" [3] "group_OP50_DAY03_vs_OP50_Day00" "group_OP50_DAY05_vs_OP50_Day00" [5] "group_OP50_DAY07_vs_OP50_Day00" "group_OP50_DAY09_vs_OP50_Day00" [7] "group_OP50_DAY11_vs_OP50_Day00" "group_SC20_Day00_vs_OP50_Day00" [9] "group_SC20_DAY01_vs_OP50_Day00" "group_SC20_DAY03_vs_OP50_Day00" [11] "group_SC20_DAY05_vs_OP50_Day00" "group_SC20_DAY07_vs_OP50_Day00" [13] "group_SC20_DAY09_vs_OP50_Day00" "group_SC20_DAY11_vs_OP50_Day00

resTC_subset_3 <- results(ddsTC1, contrast=c("group","OP50_DAY05", "SC20_DAY05"))
resTC_subset_3

alpha = 0.05
sigtab = resTC_subset_3[which(resTC_subset_3$padj < alpha), ]
sigtab
write.csv(sigtab, "OP50_Day05_vs_SC20_Day05.csv")