I have a question about the proper design for DEseq2, when there is an siControl as well as siRNA for a target gene, then there is a treated and a vehicle for both conditions.

I setup the factors and design like so:

> dds$condition <- factor(dds$condition, levels =c("siControl","siTarget"))
> dds$condition <- relevel(dds$condition, ref= "siControl") 
> dds$treatment <- factor(dds$treatment, levels = c("vehicle","treated")) 
> dds$treatment <- relevel(dds$treatment, ref ="vehicle”) 
> design(dds) <- ~ condition + treatment + condition:treatment

At first it appeared to be giving me what I thought I wanted, differences in the siTarget population upon treatment, controlling for the siControl and Vehicle condition. But then I realized that the target gene is not on the list of DE genes, even though there is a clear knockdown. I think this is because the target gene is also down in the Vehicle condition.

I can compare just siTarget vs siControl, but then this list is dominated by genes that are different in the Vehicle condition. I can also compare treated vs. vehicle, but then there is a combination of siControl and siTarget genes in the list.

Ideally what I would like to do is normalize the treated condition to the untreated condition, and then look for differences between siControl and siTarget.

Is this possible with DESeq? Is there something wrong with this approach? Is there a better way to do this analysis?

Any help or advice would be very much appreciated.

One comment: you don't need to use relevel above, if you specify the levels when you call factor. The second and fourth line aren't doing anything additional, because you already specified siControl and vehicle to be reference levels.

Can you describe what pattern you are looking for in the results? I can't tell whether you want to control for the treated vs vehicle differences or not. If you want to control for this difference, then the above design is correct, and your assay for the target gene is not working because you are not seeing a difference that is significantly different than what you see in treated vs vehicle. Am I missing something? Have you looked at plotCounts for your target?

You really need to also state how you called results, for anyone to know what results you are extracting.

You seem to be asking for two different things. If you really want genes where (control-treated/control-vehicle)/ (target-control/target-vehicle) is far from 1, then the interaction term is the way to go.

But if you want

What I would like to find is genes that are differentially expressed in the siTarget treated vs. the siControl treated.

What's simpler is to make a new column that is condition.treatment, and compare target.treated to control.treated. You do not need to make a subsetted dds of just these samples to do this. (This can be done with the interaction term in the design and some cleverness in the results call, but the condition.treatment way is much easier to follow)

If you want to compare all the targets to all the controls, while having the software understand that vehicle vs treated matters, and explains some of the variance, and should be modeled as such, you'd do ~ treatment + condition

That's why I was thinking that if the treated samples could be normalized to the vehicle samples,

This is not how DESeq works. You don't 'normalize' like that. What you can do is add treatment to the design, so it will be incorporated into the model.