Hello all,

I am new to RNA-Seq analysis and working on an experiment with three treatments compared to a control at three time points. The comparisons are structured as follows:

Reference Treatment A Treatment A+B
Vehicle day 1 vs Treatment A at day 1 vs Treatment A+B at day 1 vehicle day 3 vs Treatment A at day 3 vs Treatment A+B at day 3 Vehicle day 7 vs Treatment A at day 7 vs Treatment A+B at day 7

I created a sampleTable and design to mimic the experimental setup as follows:

sampleTable <- data.frame(
  Condition = factor(rep(c("Veh", "A", "AB"),each=9)),
  Time = factor(rep(c("1 day", "3 days", "7 days"), each = 3,time=3))
)
rownames(sampleTable) <- colnames(txi$counts)

design <- ~Condition + Time + Condition:Time

dds.temp <- DESeqDataSetFromTximport(txi,colData = sampleTable,design = design)

dds.temp2 <- dds.temp

dds.temp2$Condition <- relevel(dds.temp2$Condition,ref="Veh")

dds.temp2 <- DESeq(dds.temp2,test = "LRT",reduced = ~ Condition + Time)

data.frame(resultsNames(dds.temp2))

However, the resulting resultsNames are as follows:

resultsNames.dds.temp2.

Intercept               
Condition_AB_vs_Veh             
Condition_A_vs_Veh              
Time_3.days_vs_1.day                
Time_7.days_vs_1.day                
ConditionAB.Time3.days              
ConditionA.Time3.days               
ConditionAB.Time7.days              
ConditionA.Time7.days

I believe I am missing something, as I would like to show expression changes for each treatment compared to the control (Vehicle) at each time point. Any guidance on how to properly structure the design or extract these comparisons would be greatly appreciated!

Trying to use a conventional factor analysis, particularly if you are new, is not the way to go (I have been doing this for over 20 years now, and I don't even use a conventional factor analysis). The easier approach is to fit a cell means model and then make the contrasts you care about. To do so, you just generate a factor level for each combination of treatment and time.

> grp <- factor(gsub("\\s+", "_", apply(sampleTable, 1, paste, collapse = "_"), perl = TRUE))
> grp
 [1] Veh_1_day  Veh_1_day  Veh_1_day 
 [4] Veh_3_days Veh_3_days Veh_3_days
 [7] Veh_7_days Veh_7_days Veh_7_days
[10] A_1_day    A_1_day    A_1_day   
[13] A_3_days   A_3_days   A_3_days  
[16] A_7_days   A_7_days   A_7_days  
[19] AB_1_day   AB_1_day   AB_1_day  
[22] AB_3_days  AB_3_days  AB_3_days 
[25] AB_7_days  AB_7_days  AB_7_days 
9 Levels: A_1_day ... Veh_7_days

And now when you make your DESeq object, just use a design of ~ grp, and then you can specify whatever comparison you care about by using the contrasts argument to results. You wouldn't use an LRT test in that case though, use a Wald instead.

results(dds, c("grp", "A_1_day","Veh_1_day"))

And if you want an interaction, it's simple enough, although I have to confess I don't exactly grok the help page for contrasts and would just resort to providing a numeric vector. Let's say you want the interaction between A and AB between day 1 and 3. It's just

## check factor levels
> levels(grp)
[1] "A_1_day"    "A_3_days"  
[3] "A_7_days"   "AB_1_day"  
[5] "AB_3_days"  "AB_7_days" 
[7] "Veh_1_day"  "Veh_3_days"
[9] "Veh_7_days"
## we want (AB_3_days - AB_1_day) - (A_3_days - A_1_day)
## the coefficients for that are 1, -1, -1, 1 because without the parentheses
## it's AB_3_days - AB_1_day - A_3_days + A_1_day
## to further explain, consider this
> cbind(levels(grp), c(1,-1,0,-1,1,0,0,0,0))
      [,1]         [,2]
 [1,] "A_1_day"    "1" 
 [2,] "A_3_days"   "-1"
 [3,] "A_7_days"   "0" 
 [4,] "AB_1_day"   "-1"
 [5,] "AB_3_days"  "1" 
 [6,] "AB_7_days"  "0" 
 [7,] "Veh_1_day"  "0" 
 [8,] "Veh_3_days" "0" 
 [9,] "Veh_7_days" "0" 
## so we use
results(dds, c(1,-1,0,-1,1,0,0,0,0))

Although I think you can instead do

results(dds, list(c("AB_3_days","AB_1_day"), c("A_3_days","A_1_day")))