Deseq2 multifactor analysis

[genomica]

Hello,

I have a RNAseq experiment where I have 2 conditions and 3 replicates each and replicates are taken from different individuals. So my design is as follows:

Sample Condition Individual
S1 treated A
S2 untreated A
S3 treated B
S4 untreated B
S5 treated C
S6 untreated C

Here is my code to analyse the data using deseq2:

library('DESeq2')
directory<-"counts"
sampleFilesN1_N2 <- grep("*N[1-2]",list.files(directory),value=TRUE)
sampleConditionN1_N2<-c("treated","untreated","treated","untreated","treated","untreated")
sampleIndividual<-c("A","A","B","B","C","C")
sampleTable1<-data.frame(sampleName=sampleFilesN1_N2, fileName=sampleFilesN1_N2, condition=sampleConditionN1_N2, individual=sampleIndividual)
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable1, directory = directory,design= ~ condition + individual)
ddsHTSeq$condition <- relevel(ddsHTSeq$condition, "untreated")
dds<-DESeq(ddsHTSeq)
res<-results(dds)
res<-res[order(res$padj),]
write.csv(as.data.frame(res),file="DeSEQ2_treated_vs_untreated.csv")

However, when I run these and look at results, following message:

"log2 fold change (MAP): individual C vs A
Wald test p-value: individual C vs A"

So, its only giving me results after comparing C vs A. Am I missing anything here?

Thanks for any help.

[dpryan]

results() will just give you one of the possible comparisons (the last one in the model matrix, by default). See help(results) for how to get the results you want.

[genomica]

results() will just give you one of the possible comparisons (the last one in the model matrix, by default). See help(results) for how to get the results you want.

Thanks for the hint, that's really helpful. So, I understand if I want to get comparison of conditions considering the pairing of samples into account, I should give

res <- results(dds, contrast=c("condition","treated","untreated"))

Is this correct?

Thanks again.

[dpryan]

Yup and I see you figured out the error in your earlier attempt

[genomica]

Thanks very much Devon. Yes, I did figure out the first error

I am a bit paranoid as I am doing this analysis for the first time and hence wondered if you could confirm if my experimental design is ok.

ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable1, directory = directory,[B]design= ~ condition + individual[/B])

or do I need:

ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable1, directory = directory,[B]design= ~ condition + individual +condition:individual[/B])

Many thanks again.

[dpryan]

Go for the first one, the second one ends up directly comparing individual samples and therefore isn't reliable.

[genomica]

Thanks a lot for your help Devon.