Hello, 

First I'd like to say that I've been reading several comments and and answers from the DEVs of DESeq and I am really impressed by the amount of time you guys spend helping the community - It is really great! 

I'll briefly explain my experiment layout: 

I have two plant genotypes - V and G. At the age of 8 days they are both susceptible to a disease. At the age of 16 days only the V genotype becomes resistant while G remains susceptible.  I set G as the base genotype and 8 days as the base age - semi-arbitrarily.  

I am interested in finding genes that are deferentially expressed in V@16days compared to G@16days  but that are either:

1)  Also different in V@16days vs V8@8days - these would be genes that have changed over age in V but at 16 then are different from the G genotype.

OR 2) are also not different between the two genotypes at the 8 day susceptible age. - these would be genes that are different between genotypes specifically and only at 16 days.

I hope that was clear. I'm also not sure if those comparisons will yield the same results.

My current code:

> dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ genotype + age + genotype:age)

> dds$genotype <- relevel(dds$genotype, "G")
> dds$age <- relevel(dds$age, "8dpp")
> dds <- DESeq(dds)

> as.data.frame(colData(dds))
         genotype   age sizeFactor
Gy8-1           G  8dpp  1.5476808
Gy8-2           G  8dpp  1.2473350
Gy8-3           G  8dpp  1.0026292
Gy16-1          G 16dpp  0.8744816
Gy16-2          G 16dpp  0.9254212
Gy16-3          G 16dpp  0.9347259
Vlpk8-1         V  8dpp  0.8519697
Vlpk8-2         V  8dpp  1.1328609
Vlpk8-3         V  8dpp  0.8136137
Vlpk16-1        V 16dpp  0.7955268
Vlpk16-2        V 16dpp  0.9571078
Vlpk16-3        V 16dpp  1.2402769
> resultsNames(dds)
[1] "Intercept"          "genotype_V_vs_G"    "age_16dpp_vs_8dpp"  "genotypeV.age16dpp"

resV16G16 <- results(dds, contrast=list(c("genotype_V_vs_G","genotypeV.age16dpp")))

If I understand correctly - this should give me the results for the genotype effect at the 16day age. So basically V16 vs G16. 

How would I then modify the results to exclude or include the other contrasts?

Thanks in advance,

Ben 

PS I think this is similar to the point Simon was making in this post C: DESeq2: making specific comparisons between samples that have two conditions and from about a month ago

Hi Ben,

Excluding the sets of genes which have a small adjusted p-value for other contrasts cannot necessarily be performed using the contrast argument, but you can use setdiff() or intersect() to build the sets you want. You can build the three results tables, for the three comparisons you describe, and then filter each based on a threshold on adjusted p-value:

res1filt = subset(res, padj < 0.1)

Then setdiff(rownames(res1filt), rownames(res2filt)) will give you those genes which were significant in the first table but not in the second table. Or intersect(rownames(res1filt), rownames(res2filt)) would give you those genes which were significant in both tables. Note that lack of a small (adjusted) p-value could also be due to insufficient power, so you can't really say much about the lack of a small p-value.

Alternatively, you could test for genes which have a log fold change < theta, for some positive theta, using the alternativeHypothesis="less" argument to results. Then you would use intersection of the two tables to find genes which were DE in the first table and had a small log fold change in the second table, for a user defined value of "small". See the help in ?results and the examples in the vignette of threshold on log fold change. Note that you would have to fit the model without a prior -- DESeq(dds, betaPrior=FALSE) -- to perform this kind of significance testing looking for small log fold changes.

The 16 vs 8 day effect for the V genotype is:

results(dds, contrast=list(c("age_16dpp_vs_8dpp","genotypeV.age16dpp")))

The V vs G effect at day 8 is:

results(dds, contrast=c("genotype","V","G"))

Michael,

Thanks for your help! You are really doing a great service to the scientific community by being so responsive and helpful to us just starting to understand bioinformatics. 

Ben