This is a test version of Biostars. For the public version, visit https://www.biostars.org.
multiple comparisons in DESeq2

Hi guys,

I have the following sample information table I'm using in DESeq2:

                                      FileName        Sample Condition       Day
1           WT_2_Aligned.sortedByCoord.out.bam          WT_2        WT Uninduced
2           WT_3_Aligned.sortedByCoord.out.bam          WT_3        WT Uninduced
3           WT_4_Aligned.sortedByCoord.out.bam          WT_4        WT Uninduced
4  AAV_Con2_Day2_Aligned.sortedByCoord.out.bam AAV_Con2_Day2  Day2_Con      Day2
5  AAV_Con3_Day2_Aligned.sortedByCoord.out.bam AAV_Con3_Day2  Day2_Con      Day2
6  AAV_Con4_Day2_Aligned.sortedByCoord.out.bam AAV_Con4_Day2  Day2_Con      Day2
7  AAV_Cre1_Day2_Aligned.sortedByCoord.out.bam AAV_Cre1_Day2  Day2_Cre      Day2
8  AAV_Cre3_Day2_Aligned.sortedByCoord.out.bam AAV_Cre3_Day2  Day2_Cre      Day2
9  AAV_Cre4_Day2_Aligned.sortedByCoord.out.bam AAV_Cre4_Day2  Day2_Cre      Day2
10 AAV_Con1_Day4_Aligned.sortedByCoord.out.bam AAV_Con1_Day4  Day4_Con      Day4
11 AAV_Con2_Day4_Aligned.sortedByCoord.out.bam AAV_Con2_Day4  Day4_Con      Day4
12 AAV_Con3_Day4_Aligned.sortedByCoord.out.bam AAV_Con3_Day4  Day4_Con      Day4
13 AAV_Cre2_Day4_Aligned.sortedByCoord.out.bam AAV_Cre2_Day4  Day4_Cre      Day4
14 AAV_Cre3_Day4_Aligned.sortedByCoord.out.bam AAV_Cre3_Day4  Day4_Cre      Day4
15 AAV_Cre4_Day4_Aligned.sortedByCoord.out.bam AAV_Cre4_Day4  Day4_Cre      Day4

and I've used the following code to do my differential expression:

all(formatted_sampleinformation$Sample == colnames(formatted_countdata))

design <- as.formula(~ Condition)

model_matrix <- model.matrix(design, data = formatted_sampleinformation) 

ddsObj <- DESeqDataSetFromMatrix(countData = formatted_countdata,
                                     colData = formatted_sampleinformation,
                                     design = design)

ddsObj <- DESeq(ddsObj)

Day4_CONvsCRE <- results(ddsObj_results, contrast = c("Condition", "Day4_Cre", "Day4_Con"), alpha = 0.05)

Day2_CONvsCRE <- results(ddsObj_results, contrast = c("Condition", "Day2_Cre", "Day2_Con"), alpha = 0.05)

WTvsDay4Cre <- results(ddsObj_results, contrast = c("Condition", "Day4_Cre", "WT"), alpha = 0.05)

WTvsDay2Cre <- results(ddsObj_results, contrast = c("Condition", "Day2_Cre", "WT"), alpha = 0.05)

However, whats confusing me is that when I do:

resultsNames(ddsObj)

I get the following:

[1] "Intercept" "Condition_Day2_Cre_vs_Day2_Con" "Condition_Day4_Con_vs_Day2_Con" "Condition_Day4_Cre_vs_Day2_Con" [5] "Condition_WT_vs_Day2_Con"

leading me to believe that something may have gone wrong. Ultimately, I'd like to do the following comparisons:

WT vs Day2_Cre, WT vs Day4_Cre, Day2_Con vs Day2_Cre, Day4_Con vs Day4_Cre, Day2_Cre vs Day4_Cre

The 2nd part of the problem is that I'd like to use DESeq2's lfcShrink feature and from what I understand, you need a coef, and this can be found by running resultsNames(ddsObj), which doesn't contain some of the conditions I'd like to test (namely Day4_Con vs Day4_Cre). Have I designed this test correctly? I have a feeling I may be overcomplicating my code and have a feeling there may be a much easier way of extrapolating the comparisons I want to do.

Thanks !

deseq2 rna-seq

1 answer

If your contrasts work, you don't have to worry about what ResultsNames shows you. And you can submit a contrast to lfcshrink, see the example at the bottom

https://rdrr.io/bioc/DESeq2/man/lfcShrink.html

I don't think you need that model.matrix line in there, and I'm not sure what that all() statement is doing either.

Your four results calls need ddsObj, not ddsObj_results; the latter doesn't exist.

Log in to answer this question.