Dear Andrew,

Thanks for your detailed reply. I have already aligned all my samples using STAR and prepared count matrices with HTseq. I had performed differential expression analyses with DESeq2 separately for each experiment and looked at the intersection of the results. What I did not like of this approach is that since there are two experiments and 3 treatment groups each, I end up with so many pairwise comparisons that are messy to summarize, and also the two experiments for some reason seemed to have different detection power overall so it's complicated to assess the overlap of the same pairwise comparisons across experiments just based on a cutoff of adjusted P values, and this applies also for downstream analyses such as GO enrichment. This is why I was looking into a way to analyze the two experiments in one go (ending up with less comparisons) and see how the results would look like compared to the previous approach. But if I cannot include my two control groups then it's useless (as they are the ones I'm especially interested in comparing the treated samples to). Everything in the experiments is the same (experimental procedures and sequencing approach) so that should not be a problem, but indeed PCA analyses show there is a batch effect of the two experiments. The only thing that differed was the way control samples were treated. What I'm now thinking (if really there is no way of including control groups in an overall analysis) is to merge the controls and treat them as they were the same (conceptually it wouldn't be too far-fetched as really I want to know what in the treatments is different from any of those controls, but I'm not sure if I'd be violating any important assumption). What would you think about that?

Thanks again for your help!