Hi all,

I study the effect of cadmium (= a heavy metal) on leaf growth.

To give you some background on the experiment:
To study the leaf in greater detail, three distinct zones of the leaf are selected: mer, elo and mat (each known for a specific biological process, without bothering now about the details). 
I work with three treatments: control, mild and severe (where control = no cadmium, mild = low cadmium dose, severe = high cadmium dose)
For this experiment, I’ve performed a transcriptome study with three biological replicates for each possible combination (e.g. 3xcontrol mer, 3xcontrol elo, 3xcontrol mat, 3xmild mer, etc.)

Experimental setup in short:
Three cadmium treatments: control, mild, severe
Three zones in the leaf: mer, elo and mat

I’ve already made some progress:
- Count tables (unique gene reads) made using CLC
- Read in the count table in R and got it running in the DESeq2 package
- Made pair-wise comparisons using the group-functionality in DESeq2:

# construction to compare specific groups which are combinations of variables
dds1$group <- factor(paste0(dds1$treatment, dds1$zone))
design(dds1) <- ~ group
dds1 <- DESeq(dds1)
resultsNames(dds1)
resu1 <- results(dds1, contrast=c("group", "controlmer", "severemer"))
head(resu1)
write.table(resu1, "/Users/…/controlmer_vs_severemer.txt", sep="\t")

Where I’m stuck:
From past statistics courses (on the simple two-way ANOVA), I’m used to obtaining ONE p-value for treatment effect, ONE p-vale for zone effect and ONE p-value for the interaction between treatment and zone. I’m also interested in the effect of my treatment, zone and the interaction for each gene in my transcriptome study, but I can’t figure out how to obtain this ONE p-value for each (treatment, zone and the interaction).

My code in R:

dds1 <- DESeqDataSetFromMatrix(countData = mydata, colData = colData, design = ~ zone + treatment + zone:treatment)
dds1 <- DESeq(dds1)
resultsNames(dds1)
res1 <- results(dds1)
head(res1)

Output of “resultsNames(dds1)”:

[1] "Intercept"    "zone_mat_vs_elo"   "zone_mer_vs_elo"   "treatment_mild_vs_control"   "treatment_severe_vs_control"   "zonemat.treatmentmild"   "zonemer.treatmentmild"       "zonemat.treatmentsevere"   "zonemer.treatmentsevere"    

I only manage to obtain p-values for specific contrasts, for instance:

# treatment mild vs control
resultsNames(dds1)
results(dds1, contrast=c("treatment", "mild","control"))
results(dds1, list( c("treatment_mild_vs_control") ))

# treatment severe vs control
resultsNames(dds1)
results(dds1, contrast=c("treatment", "severe","control"))
results(dds1, list( c("treatment_severe_vs_control") ))

Question:
How can I obtain a p-value for the overall treatment effect, a p-value for the overall zone effect and a p-value for the interaction between treatment and zone for each gene individually using DESeq2? 

I think that, obtaining these p-values, will direct me to the genes of interest (especially the interaction, where the zone effect is not consistent across treatments). (Feel free to correct me when I'm wrong by putting to much emphasis on these p-values. Maybe I should only focus on the contrasts (using the group-functionality in DESeq2 as shown above).)

Thank you in advance!

Kind regards,
Jonas

University of Antwerp
Belgium

hi Jonas,

When you run anova() in R, the function builds a series of nested models, by removing terms one-by-one.

> y <- rnorm(40)
> x <- factor(rep(1:2,each=20))
> z <- factor(rep(c(1,2,1,2),each=10))
> anova(lm(y ~ x + z + x:z))
Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value Pr(>F)
x          1  0.7034 0.70339  0.8637 0.3589
z          1  0.0838 0.08376  0.1028 0.7503
x:z        1  1.9524 1.95238  2.3972 0.1303
Residuals 36 29.3195 0.81443

You can explicitly specify the models and get the same sum of square terms

> anova(lm(y ~ 1), lm(y ~ x), lm(y ~ x + z), lm(y ~ x + z + x:z))
Analysis of Variance Table

Model 1: y ~ 1
Model 2: y ~ x
Model 3: y ~ x + z
Model 4: y ~ x + z + x:z
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     39 32.059
2     38 31.356  1   0.70339 0.8637 0.3589
3     37 31.272  1   0.08376 0.1028 0.7503
4     36 29.319  1   1.95238 2.3972 0.1303

You can do something similar with DESeq2, by manually specifying the two models to compare and perform a likelihood ratio test, see the vignette section on the LRT.

I will note that anova() does something slightly different if you specify only two models as opposed to the full sequence (other than the last pair) because it is comparing to the residuals from the largest model provided.

> anova(lm(y ~ 1), lm(y ~ x))
Analysis of Variance Table

Model 1: y ~ 1
Model 2: y ~ x
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     39 32.059
2     38 31.356  1   0.70339 0.8524 0.3617