Dear DESeq2 developers & Bioconductor community.

I was hoping somebody could clarify/confirm on how to use the following experiment data in DESeq2:

1. I have pairwise expression data from two distinct loci, A and B (the position), and under two different conditions (treated and untreated).
2. I am interested in exploring whether there is a difference in the ratio A/B between the two conditions. Just to clarify: I am not interested in finding out whether there are significant differences between A and B (there are), but whether those differences between A and B change between condition treated vs. untreated.

The data is available in quadruplicates(untreated)/duplicates(treated):

> head(countData)
    A_treated_rep1 A_treated_rep2 A_untreated_rep1 A_untreated_rep2
1            48             3             29             19
2           123            35             34             50
3            96            36             18             44
4            61            39             25             14
5           861           715            432            249
6            99            48             93             20

  A_untreated_rep3 A_untreated_rep4 B_treated_rep1 B_treated_rep2
1             24             25            49            30
2             11             61            75            14
3             22             51            24            17
4             14             28            56            11
5            134            118           757           564
6             12             67            56            11

  B_untreated_rep1 B_untreated_rep2 B_untreated_rep3 B_untreated_rep4
1              0              0              0              0
2              0              0              0              0
3              0              0              0              0
4              0              0              0              0
5             65             47             22             29
6              0              0              0              0

> colData
                  position condition
A_treated_rep1       A      treated
A_treated_rep2       A      treated
A_untreated_rep1     A      untreated
A_untreated_rep2     A      untreated
A_untreated_rep3     A      untreated
A_untreated_rep4     A      untreated
B_treated_rep1       B      treated
B_treated_rep2       B      treated
B_untreated_rep1     B      untreated
B_untreated_rep2     B      untreated
B_untreated_rep3     B      untreated
B_untreated_rep4     B      untreated

I proceed in the following way:

dds<-DESeqDataSetFromMatrix(countData=countData,colData=colData,design = ~condition+position+condition:position);
dds<-DESeq(dds);

Now comes the part, where I am not so sure as to what I am doing is correct; I look at the contrast

res<-results(dds,name="conditionuntreated.positionB");

As far as I gather from the manual, different forum posts here and the ?results help page, this should tell me if there is difference in the overall position effect (B/A) between condition treated and untreated.

Is this correct?

Furthermore, how would I go and interpret the log2 FC in this case? Does the log2 FC directly translate into the difference log2(B/A(treated)) - log2(B/A(untreated))?

And one more thing: Would a design = ~condition:position work as well, i.e. keeping only the interaction term in the glm? I guess this would have the effects due to condition and position be absorbed in the coefficient for the interaction term (at least to some degree). Which probably would not be too sensible?

I would appreciate any guidance/insight.

"Furthermore, how would I go and interpret the log2 FC in this case? Does the log2 FC directly translate into the difference log2(B/A(treated)) - log2(B/A(untreated))?"

Yes. Except that the sign is opposite: log2(B/A(untreated)) - log2(B/A(treated))

This is because the result name, "conditionuntreated.positionB", contains the factor levels that are in the numerator, i.e., counted positive (and, by default, R orders factor levels alphabetically, placing "untreated" after "treated").