Hi,

I am using DESeq2 to analyze differential expressed genes (DEGs) from counts generated by Illumina sequencing.

I have a 3-factor experiment:

• genotype (with 3 levels): A, B, and C, each with 6 biological replicates (3 for each treatment)
• treatment (with 2 level): mDr (the treatment) and WW (the control)
• timepoint (with 3 levels): T1, T2 and T3.

In order to analyze genes that are influenced by an individual component and those influenced by two or all variables, I'm interested in the following results/questions:

1. How many DEGs are influenced by an individual component (genotype, treatment, timepoint)?
2. How many DEGs are influenced by the combination of two or all components (genotype:treatment + genotype:timepoint + treatment:timepoint + genotype:treatment:timepoint)?.

So, I tried as follows:

library(DESeq2)

# Load row counts
countData<-read.delim("all_counts_matrix_T1_T2_T3", sep="\t", header=TRUE, row.names = 1)
col_ordering = c(1,3,5,2,4,6,7,9,11,8,10,12,13,15,17,14,16,18, 19,21,23,20,22,24,25,27,29,26,28,30,31,33,35,32,34,36,37,39,41,38,40,42,43,45,47,44,46,48,49,51,53,50,52,54)
countData_ord = countData[,col_ordering]

# Create countData as matrix
countData_m <- as.matrix(countData_ord)
# Interger numbers in matrix
storage.mode(countData_m) = "integer"

#Create sample information
colData<-read.delim("colData.txt", sep="\t", header=TRUE, row.names = 1)

head(countData_m)
colData

rownames(colData) == colnames(countData_m)

# Count matrix and sample information input and design formula
dds<-DESeqDataSetFromMatrix(countData = countData_m, colData = colData, design = ~ genotype + treatment + timepoint + genotype:treatment + genotype:timepoint + treatment:timepoint + genotype:treatment:timepoint)

dds$treatment<- relevel(dds$treatment, "WW")

dds$treatment
dds$genotype
dds$timepoint

## Filtering to remove rows with 0 reads or only a single count across all samples
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)
dds

# Run DESeq2
dds = DESeq(dds)
# Calling results
res <- results(dds)
mcols(res, use.names=TRUE)

# Print the coefficients
resultsNames(dds)

summary(res)
write.table(res, file='genotype_treatment_timepoint/res_04082016.txt', col.names=T, row.names=T, quote=F, sep="\t")

In resultsNames(dds) I get:
[1] "Intercept"                          "genotype_B_vs_A"                    "genotype_C_vs_A"                    "treatment_mDr_vs_WW"             
[5] "timepoint_T2_vs_T1"                 "timepoint_T3_vs_T1"                 "genotypeB.treatmentmDr"             "genotypeC.treatmentmDr"          
[9] "genotypeB.timepointT2"              "genotypeC.timepointT2"              "genotypeB.timepointT3"              "genotypeC.timepointT3"           
[13] "treatmentmDr.timepointT2"           "treatmentmDr.timepointT3"           "genotypeB.treatmentmDr.timepointT2" "genotypeC.treatmentmDr.timepointT2"
[17] "genotypeB.treatmentmDr.timepointT3" "genotypeC.treatmentmDr.timepointT3"

First of all, this is correct?

And then, how can I extract the specific results I wanted (question 1 and 2)? I've tried to read up on this, but haven't been able to find  how to code this.

Any guidance is much appreciated. 

Thanks in advance

Chiara

If you really want to test inclusion of specific factors, rather than looking at the contrast between levels within a factor, then I'd suggest using a likelihood ratio test, which DESeq2 conveniently provides:

results(dds, type="LRT", full= ~genotype + treatment + timepoint, reduced= ~genotype + treatment)

will inform you which genes need to take 'treatment' into account in order to explain the variability between samples, over and above the variability induced by genotype and treatment.  You may want to consider the full=~treatment, reduced=~1 comparison to see if the null model you choose alters the output.  You can do this for the other factors (include and exclude it from the full and reduced models) to examine their influence, and similarly you can add second (and third) order interactions:

results(dds, type="LRT", full= ~genotype + treatment + timepoint + genotype:treatment, reduced= ~genotype + treatment + timepoint)

to look at question 2 (in this case looking at genotype:treatment interactions, but the others are similar).  For the timepoint comparison, it might also be worth adding a 'biological replicate' factor, so you could see if individual replicates are behaving consistently. 

This is a complex experimental design and there are many ways to pose hypotheses here and to test them with DESeq2. I'd recommend  you partner with a local statistician to help pose specific hypotheses as contrasts or likelihood ratio tests and to help interpret results. There is nothing specific about DESeq2 for the comparisons you'd like to make, you would use the same design and contrasts as with a normal linear model.