DESeq2 on multiple comparison.

Question

[a.archana]

I am analysing data generated on Arabidopsis plants. My experiment design is as follows:

  condition treatment
1 FDR       Female_insect
2 SAX       Female_insect
3 FDR       Male_insect
4 SAX       Male_insect
5 FDR       No_insect
6 SAX       No_insect

The question I am trying to answer is: There are some genes which are known to induce in SAX plant compared to FDR in male treatment, but not induced in SAX compared to FDR in No insect (control) and female treatment.

I understand I can use following command in DESeq2

dds <- DESeqDataSetFromMatrix(countData = dataCountTable, colData = dataDesign, design = ~ condition)
dds <- estimateSizeFactors(dds)
countTable <- counts(dds, normalized=TRUE) 
dds <-DESeq(dds) 
dds$group <- factor(paste0(dds$treatment, dds$condition))
design(dds) <- ~ group
dds <- DESeq(dds)
resultsNames(dds)
results(dds, contrast=c("group", "SAX", "male_insect"))

Similarly I can get DE genes for all other comparison. And can do venn diagram for overlap?

Is there any other way to answer this question. And can anyone please comment on the command above if this is right?

Any help would be appreciated. Thanks in advance.

Answer 1

[Devon Ryan]

I presume you have more replicates than this, otherwise you have no reliable way of answering the question.

I think what you really want is a mode like condition*treatment and to look at differentially expressed genes (presumably with a positive fold-change) in the Male_insect:SAX interaction. That's a bit easier to specify with a factorial design than trying to use a group design and using contrasts.

[a.archana]

Thanks Devon. I got 4 replicated for all of the above combinations.

Could you please comment a bit more on factorial design? I am not very sure what you mean. Thanks

[Devon Ryan]

Instead of DESeqDataSetFromMatrix(countData = dataCountTable, colData = dataDesign, design = ~ condition), try DESeqDataSetFromMatrix(countData = dataCountTable, colData = dataDesign, design = ~ treatment*condition) instead. This is then a model with both condition and treatment as independent main effects as well as their interactions. In many cases this sort of design is a bit easier to deal with, especially when what you want is to look at an interaction between main effects.

[a.archana]

Thanks Devon. Got it.

[a.archana]

Dear Ryan,

Thanks for your time and patience.

Following your comment I added following to my design:

dds <- DESeqDataSetFromMatrix(countData = dataCountTable, colData = dataDesign, design = ~ genotype+treatment+genotype:treatment)

dataDesign genotype treatment sample1 F FDR sample2 F FDR sample3 F FDR sample4 F FDR sample5 F SAX sample6 F SAX sample7 F SAX sample8 F SAX sample9 M FDR sample10 M FDR sample11 M FDR sample12 M FDR sample13 M SAX sample14 M SAX sample15 M SAX sample16 M SAX sample17 C FDR sample18 C FDR sample19 C FDR sample20 C FDR sample21 C SAX sample22 C SAX sample23 C SAX

dds <- DESeq(dds)

resultsNames(dds)

[1] "Intercept" "genotype_F_vs_C" "genotype_M_vs_C"
[4] "treatment_SAX_vs_FDR" "genotypeF.treatmentSAX" "genotypeM.treatmentSAX"

1. I do not understand : genotypeF.treatmentSAX and genotypeM.treatmentSAX interaction meaning?

2: How can I answer this question: Which genes are DE in male-exposed SAX versus male-exposed FDR plants, that are not DE in other comparisons. Meaning these genes should be not DE or significantly less DE or DE in the opposite direction in female-exposed FDR versus female-exposed SAX plants, or in male-exposed FDR or control FDR plants?

I hope its not too much of your time. many thanks

[Devon Ryan]

For question 2 it's easier to use a non-factorial model like you showed in your original post. For the interaction, let's take the example of genotypeM.treatmentSAX. These will be the genes that change more than one would expect due to a combination of an M genotype or the SAX treatment together. That is, they change more than expected if the genotype and treatment act independently, which suggests that there's a biological level at which they're interacting with each other. The up-regulated genes in this are going to be the ones you're most interested in (though you might need to ensure that they're not also up-regulated in genotypeF.treatmentSAX).

[a.archana]

Thanks Devon. Interaction is clear to me now :)

I will try the first question and get back to you. Thanks for your time! thanks a lot.