I am working on analyzing a dataset using DESeq2. The main goal of the analysis is to compare the gene expression profiles between treated and untreated group. Since we also have sex information on the patient, we are also looking at sex-specific differences of the main condition(Treatment). I have added a code snippet below where I was able to compare the treatment effects in males and females separately. I have come across new information on the patients about their heart disease status.
Now my question is what is the best way to write deseq2 design and contrast in order to extract something like, what are the differences between treatment groups for men who also have a heart disease.
PatientID Condition Sex Heart disease ABC1234 Treated M Yes ABC1235 Treated F No ABC1236 Treated M Yes ABC1237 Treated F Yes ABC1238 Treated M Yes ABC1239 Untreated F No ABC1240 Untreated F Yes ABC1241 Untreated M No ABC1242 Untreated F Yes Gr <- factor(paste0(Pheno$Condition,Pheno$Sex)) dds <- DESeqDataSetFromMatrix(countData = myFile, colData = Pheno, design= ~Gr) dds <- DESeq(dds,parallel=TRUE) resultsNames(dds)  "Intercept"  "Gr_UntreatedM_vs_UntreatedF"  "Gr_TreatedF_vs_UntreatedF"  "Gr_TreatedM_vs_UntreatedF" Male <- results(dds, contrast=c('Gr','TreatedM','UntreatedM')) Female <- results(dds, contrast=c('Gr','TreatedF','UntreatedF'))