Hi,

I am doing DE analysis on three breed groups (FS, TX, F1) with two diet conditions (C, F, where C is the normal diet and F is the additional diet). I want to find DE genes between breeds (FS vs TX, FS vs F1 and TX vs F1) both with diet effect and without diet effect. In addition, I want to find that are differentially expressed within breed due to diet (FSc vs FSf, TXc vs TXf, F1c vs F1f). I designed the experiment like this:

sampleFiles<-list.files(directory, pattern = "*_counts")
sampleinfo<-read.csv("SampleTable_mRNA_all.csv", header=TRUE)
row.names(sampleinfo)
sample<-sampleinfo$Sample
breed<-sampleinfo$Breed
diet<-sampleinfo$Diet
id<-sampleinfo$SampleID
sampletable<-data.frame(sampleName=sample, fileName=sampleFiles,breed=breed, diet=diet, id=id)
group<-factor(paste(breed, diet, sep="."))
sampletable<-cbind(sampletable, group=group)
ddsMat<-DESeqDataSetFromHTSeqCount(sampleTable = sampletable, directory = directory, design=~group)
ddsMatcol<-collapseReplicates(ddsMat, sampletable$id)
dds<-DESeq(ddsMatcol)
res<-results(dds)

Then I made comparisions like this:

For example, within-breed diet effect on TX:
dTX_Diet<-results(dds, contrast=c("group", "TX.F", "TX.C"))
dTX_Diet_sig<-subset(dTX_Diet, padj< 0.05)

similarly, for identifying DE genes between FS and TX due to diet effect

dFS_TX<-results(dds, contrast=list(c("groupFS.F", "groupTX.C"), c("groupFS.C", "groupTX.F")))
dFS_TX_sig<-subset(dFS_TX, padj<0.05)

Could anybody please tell me if I would rather need to use LRT method?

In addition, I am interested to find DE genes in FS vs TX, FS vs F1 and TX vs F1 without including diet factor, in this case, how can I exclude the diet effect or the only way is to compare is that between FSc vs TXc and so on? 

Kisun

hi Kisun,

The way you have set up the design is good for comparing pairs of groups, for example, comparing within breed F vs C is correct as you have it.

Also the difference of differences contrast looks correct to me, testing the diet effect across groups, (FS.F - FS.C) - (TX.F - TX.C), by rearranging to: FS.F + TX.C - (FS.C + TX.F).

In DESeq2, we use Wald tests for such contrasts, so your current approach is what I would recommend.

If you wanted to find DE genes across groups controlling for diet but not specifically testing across groups within one diet, then I would run DESeq() with a separate design:

~ diet + breed

And then test on breed, using, e.g. contrast=c("breed","X","Y").

The reason I would re-run with this design, is that typically, if I read that a model "controlled for diet" I would assume that it involves less terms than your former design. I would assume that you meant that the diet effect was modeled as a single multiplicative fold change for each gene, but that it does not interact with breed group.