DESEq2 Factor/Level Design

Question

[sdbaney]

I have 3 muscle types

1. keel
2. mantle
3. nuchal

I have 3 replicates between each and want to perform differential expression across all 3 to see which transcripts are being significantly differentially expressed. I have been using Cuffdiff because it's simple and the output is nice in that it has pairwise comparisons between everything.

I tried this with DESeq2, factor being muscle, level 1 being keel (select the 3 replicate count files), level 2 mantle, level 3 nuchal.

But when I ran the program and looked at the output it reported significance and expression values but it didn't tell me what samples it belonged to. I looked more into DESeq2 and it seems that the first level is accepted as the reference and then any further levels are just compared back to that. So I would have keel as the reference and then only get output of keel vs mantle, and keel vs nuchal. However I still don't know which is what since it doesn't tell me the samples it originates from. Additionally, I don't have a mantle vs nuchal comparison.

I'm assuming there's just something I am not understanding in the designing aspect of this setup. Can anyone give me some guidance?

Answer 1

[Asaf]

You can get all the comparisons using the results() function, just put the right levels in the contrast. See http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts

Answer 2

[h.mon]

Please show your DESeq2 code. But basically you want to use the contrast argument from the results() function:

res_mantle_vs_nuchal <- results(dds, contrast=c("muscle","mantle","nuchal"))
res_keel_vs_nuchal <- results(dds, contrast=c("muscle","keel","nuchal"))

Answer 3

[blawney]

What do you mean by "I still don't know which is what since it doesn't tell me the samples it originates from?"

The table produced by the results function (as given in the other answer), produces where each row gives information about a single gene (e.g. fold change, p-val, adjusted p-val). The table tells you "on average" (colloquial, non-technical usage) if the tissues have different expression. It's not going to tell you which of the the samples are themselves individually "significant" (whatever that might mean).

If you want to examine the results (always good to spot-check, especially if your designs are complex), you might use the counts(dds, normalized=TRUE ) function to extract the table of normalized counts. Then plot those normalized expressions to see that they are in "visual agreement" with the DESeq2 findings.

[sdbaney]

What I meant was, what sample - which tissue has the higher/lower expression? It doesn't tell me for example, that keel is expressing a value of x.