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)
     [1] "Intercept"                            
     [2] "Gr_UntreatedM_vs_UntreatedF"
     [3] "Gr_TreatedF_vs_UntreatedF"  
     [4] "Gr_TreatedM_vs_UntreatedF" 

    Male <- results(dds, contrast=c('Gr','TreatedM','UntreatedM'))
    Female <- results(dds, contrast=c('Gr','TreatedF','UntreatedF'))

How many samples do you have in total? If you have the samples listed above, this is not sufficient to answer, "what are the differences between treatment groups for men who also have a heart disease". If you have many more samples, such that you have sufficient numbers in each combination, you could use the same approach as above but combining also heart disease status.