Hi everyone,

I am trying to carry out some analysis for an experiment which has three dimensions (treatment, genotype and time) and while I have managed to extract the pairwise comparisons I am interested in, I would like to go further and look at the genes which behave differently upon my drug treatment depending on either the genotype or the time factor (for a given genotype) - I believe that this is the interaction term. However, I am not entirely sure how to do this. I have read some of the posts here, the vignette and the examples in ?results, but would like to check here.

To describe my experiment in detail: I have two genotypes (MYC On and MYC Off), which were treated with vehicle or drug (Veh and Drug), and the mRNA levels measured after 2, 6 and 18 hours of exposure.

To obtain the pairwise comparisons, I set the design using a grouping variable, and then extracted each comparison of interest:

dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~Genotype + Time + Treatment)
dds$group <- factor(paste0(dds$Genotype, dds$Time, dds$Treatment))
design(dds) <- ~group

#Extract each pairwise comparison, e.g.
MYC_ON_18h_TreatedvsWT <- lfcShrink(dds, contrast=c("group", "ONT18DRUG", "ONT18VEH"), type="ashr")

This worked well, but as stated, I would now like to see which which genes are differentially affected by my treatment based on the:

a) Genotype (for each timepoint, if possible).
b) Time (for a given genotype).

I am unsure as to how to write the code to pull this information out. I am under the impression that I cannot do this anymore using a grouping variable, and have to define the design formula explicitly, based on the phrasing on ?results. I am not sure how then to extract the relevant data for each timepoint and comparison. I have written something which I think might work, but would appreciate another opinion.

dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~Genotype + Time + Treatment + Genotype:Treatment + Time:Treatment)

#check reference levels
dds$Genotype
#Answer: Levels: ON OFF 
dds$Time
#Answer: T18 T2 T6
dds$Treatment
#Answer: DMSO DRUG
#Therefore current reference level is MYC ON, vehicle, 18h.
#It would make more sense to me from an experimental flow standpoint if the current reference level is MYC ON, vehicle, 2h, so I change this below.

dds$Time = relevel(dds$Time, "T2")

dds <- DESeq(dds)

#Check the comparisons made
resultsNames(dds)

#[1] "Intercept"              "Genotype_OFF_vs_ON"          "Time_T18_vs_T2"         "Time_T6_vs_T2"         
#[5] "Treatment_DRUG_vs_DMSO" "Genotype_OFF.TreatmentDRUG"   "TimeT18.TreatmentDRUG"  "TimeT6.TreatmentDRUG"  

#Extract results for 2h timepoint - differential effect of drug based on genotype.
results2h <- results(dds, name="Genotype_OFF.TreatmentDRUG")

#Extract results for 6h timepoint - differential effect of drug based on genotype.

#Extract results for 18h timepoint - differential effect of drug based on genotype.

#Extract results for MYC on - differential effect of drug based on time.
2v6 <- results(dds, name=" TimeT6.TreatmentDRUG")  #(2h vs 6h, MYC ON)
2v18 <- results(dds, name=" TimeT18.TreatmentDRUG") #(2h vs 18h, MYC ON)
6v18 <- results(dds, contrast=list("TimeT18.TreatmentDRUG", "TimeT6.TreatmentDRUG") #(6h vs 18h, MYC ON).


#Extract results for MYC off - differential effect of drug based on time.

However, I am not sure how to extract the resutls for the 6h timepoint - differential effect of drug based on genotype, 18h timepoint - differential effect of drug based on genotype or the results for MYC off - differential effect of drug based on time. Would I have to re-define the reference level and then rerun dds? Example below:

dds$Time = relevel(dds$Time, "T6")
dds <- DESeq(dds)

#Extract results for 6h timepoint - differential effect of drug based on genotype.
results6h = results(dds, name="Genotype_OFF.TreatmentDRUG")

dds$Time = relevel(dds$Time, "T18")
dds <- DESeq(dds)

#Extract results for 18h timepoint - differential effect of drug based on genotype.
results18h = results(dds, name="Genotype_OFF.TreatmentDRUG")

Any help would be greatly appreciated.