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!
2 answers
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.
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")
Log in to answer this question.