DESeq2 extract desired result

Question

[beans]

Hi everyone, I have a bulk RNAseq experimental design with two genotypes (WT and Mutant), and two conditions (Untreated and Treated). There are two sets of genes that I hope to find:

1. Genes with expression patterns that change differently due to treatment across genotypes (eg. genes that increased expression with treatment in WT but decreased expression with treatment in Mutant)

I did LRT where the full model is Genotype + Condition + Genotype*Condition, and compared to the reduced model Genotype + Condition.

1. Genes with expression patterns that change similarly due to treatment across genotypes (eg. genes that increased expression with treatment in WT and also increased expression with treatment in Mutant)

I am a little lost on how to extract results to find this subset of genes. I have thought about merging the sets of genes that are differentially expressed in two pairwise comparisons: Genes that are differentially expressed when treated within WT samples - results(dds, contrast=c("condition","Trt","Ctrl")) AND Genes that are differentially expressed when treated within KO samples - results(dds, list( c("condition_Trt_vs_Ctrl","genotypeMU.conditionTrt") ))

Thanks!

[H Mazumder]

Do you have an update on this work? Was the provided answer useful? Thanks.

Answer 1

[H Mazumder]

Hi beans,

You may try a simpler design which also accounts for Genotype*Condition interaction. You can combine your two variables (treatment and genotype) into one variable (say, group). Then you can directly call the desired comparisons using the result() function of DESeq2. Here, is an example code that you can add in your existing DESeq object dds:

dds$group = factor(paste0(dds$genotype, dds$condition))

This gives you a group variable with labels: WTUntreated, WTTreated, MutantTreated, MutantUntreated. Next update the deseq2 object dds.

design(dds) = ~ group
dds = DESeq(dds)

#extract the results and merge into one data frame:

res1 = as.data.frame(results(dds, contrast=c("group", "WTTreated", "WTUntreated")))
res1$gene_name = rownames(res1)
colnames(res1)[1:6] = paste(colnames(res1)[1:6], "WT", sep="_")
res2 = as.data.frame(results(dds, contrast=c("group", "MutantTreated", "MutantUnreated")))
res2$gene_name = rownames(res2)
colnames(res2)[1:6] = paste(colnames(res2)[1:6], "Mutant", sep="_")

res = merge(res1, res2, by = "gene_name") #merge

#For question 1, the genes are in this data set:
res.q1 = res[res$log2FoldChange_WT>0 & res$log2FoldChange_Mutant<0, ]
res.q1.sig = res.q1[res.q1$padj_WT < 0.05 & res.q1$padj_Mutant < 0.05, ]

#For question 2, the genes are in this data set:
res.q2 = res[res$log2FoldChange_WT>0 & res$log2FoldChange_Mutant>0, ]
res.q2.sig = res.q2[res.q2$padj_WT < 0.05 & res.q2$padj_Mutant < 0.05, ] 

res.q1.sig gives genes that were upregulated in WT group but downregulated in Mutant group due to treatment. res.q2.sig will give you genes that were upregulated due to treatment in both WT and Mutant groups. The above approach is consistent with DESeq2 vignettes.

I hope this answer helps. Thanks!

~H. Mazumder