[Michael Love]

hi Gwenola,

Check in the Examples section of the help for ?results:

You should find, Example 2: two conditions, two groups, with interaction term.

[gtosser]

hi gamb,

Firstly, as you are just starting to analyze, best to update to Bioconductor 3.0 and DESeq2 version 1.6 which was just released.

As the analysis is a bit complicated, I'd recommend using DESeq(dds, betaPrior=FALSE) as this simplifies the building of contrasts with interactions.

The design should be

~ tissue*chem1*chem2

...because you're interested in both the interactions of chem1 and chem2, and in the specific effects for each tissue.

The list of effects is given by

resultsNames(dds)

You can build up the contrasts using a list like below. I'll give some examples and see if you can extend to the others. Comparing the tissues alone is simple by specifying the name of the factor and the two levels to compare, e.g.:

results(dds, contrast=c("tissue","C","B"))

Then some basics about interaction terms: the effect of chem 1 in tissue B will be:

results(dds, contrast=list(c("chem1","chem1.tissueB")))

...because you need to add the main chem 1 effect (the effect for tissue A, the base level of tissue) and the extra chem 1 effect for tissue B.

If you wanted to test if chem 1 has an effect in tissue B different than tissue A you would only test the interaction term "chem1.tissueB", because this is that difference.

Tissue A is the base level, so the effect of chem 1 in tissue A is simply

results(dds, name="chem1")

The effect of chem 1 and chem 2 in tissue B is a bit lengthy, you need to add main effects and all relevant interaction terms:

results(dds, contrast=list(c("chem1","chem2","chem1.tissueB","chem2.tissueB","chem1.chem2","chem1.chem2.tissueB")))

It's simpler for tissue A because this is the base level:

results(dds, contrast=list(c("chem1","chem2","chem1.chem2")))

And if you want to test only the interaction, answering if the effect of the two chems is different than each separately, you remove some of these effects. For tissue B this contrast would be:

results(dds, contrast=list(c("chem1.chem2","chem1.chem2.tissueB")))

For tissue A it is again simpler:

results(dds, name="chem1.chem2")

Note that another way to compare individual groups would be to create a new variable which describes all groups:

dds$group = factor(paste(dds$chem1, dds$chem2, dds$tissue))
design(dds) = ~ group

Then use the following to compare to specific groups:

results(dds, contrast=c("group","...","..."))

Dear Michael,

I am afraid I cannot build my contrsats either.
I have cells from 2 species (bovin, porc) that have been treated or not (treatment=pos or treatment=neg)

Here is my model
design (ddsX)<-~species+treatment+species:treatment

resultsNames(ddsX)
#[1] "Intercept" "species_porc_vs_bovin" #"treatment_pos_vs_neg" "speciesporc.treatmentpos"

If I understand well, I can find the genes with a treatment effect with:
results(ddsX, contrast=c("treatment","neg","pos"))

then the genes with a species effect with:
results(ddsX, contrast=c("species","bovin","porc"))

But how can I find the genes:
- from species "bovin" with a treatment effect
- from species "porc" with a treatment effect
- with interaction (that do not react to the treatment in a same way in "bovin" or "porc")?

I am using DESeq2_1.6.2

Thanks for your help

Gwenola

[Michael Love]

It might be a problem in .RData, which is loaded every time. Try deleting this file. (although doesn't explain why you get it on two separate computers)

[gamb]

Thats the thing... I am running the latest GUI

[Michael Love]

I'm not sure about this, but I would recommend following the message and updating R to the latest GUI version.

[gamb]

So I got R to throw up this error (not a segfault...these don't happen too often) right after resultsNames(dds)

2014-10-27 19:31:47.892 R[5511:707] *** RController: caught ObjC exception while processing system events. Update to the latest GUI version and consider reporting this properly (see FAQ) if it persists and is not known. 
*** reason: *** -[__NSArrayM objectAtIndex:]: index 18446744073709551615 beyond bounds [0 .. 85]
*** name: NSRangeException, info: (null)
*** Version: R 3.1.1 () R.app R 3.1.1 GUI 1.65 Mavericks build(null)
Consider saving your work soon in case this develops into a problem.

[Michael Love]

Yeah, that won't be from DESeq2 then. And maybe ask the phyloseq maintainers if they have seen this on their github site, although their data structures are data.frame and matrix so I wouldn't know why there would be a problem there. If you are loading objects from memory at the beginning those could be corrupted.

[gamb]

What do you mean "at will"? After running a particular function? At what points in the analysis to you get these?

It's been pretty random when they happen- usually when I'm not running code (often just editing script files prior to running them) so I haven't been able to track down what particular thing is triggering them. All I do know is that they happen in sessions when I've loaded those packages. I think it might be tied into the phyloseq data structures.

[Michael Love]

A safer way to reorder columns would be to build the DESeqDataSet with a design of ~ 1, then rearrange the columns of the DESeqDataSet, then add a new column anim.nested, then add the desired design. What do you mean "at will"? After running a particular function? At what points in the analysis to you get these?

[gamb]

Dang it! Good catch, my columns did not match up.

When I get another error I'll screenshot it and post it here. They seem to pop-up at will. I have a normal release of R and the errors occur at my lab computer + personal computer.

[Michael Love]

The geoMeans code you have is the old one, where I forgot to exponentiate the mean of log counts. See the link in my previous reply.

Otherwise seems right. Be very careful that the columns of your count matrix line up with colData.

I've never seen or heard of a segfault from our code. Do you have a normal release of R or did you build yourself?

[gamb]

Ok! So, just checking in to make sure this is right... I reformatted the data as described named "d", "bbb" is a phyloseq data object, and "m" is the matrix model.

counts=as.data.frame(otu_table(bbb))
dds=DESeqDataSetFromMatrix(countData=counts, colData=d, design=~ chem1 + chem2 + chem1:chem2 + chem1:chem2:anim.nested + chem1:chem2:tiss)

cts = counts(dds)
geoMeans = apply(cts, 1, function(row) mean(log(row[row != 0])))
dds = estimateSizeFactors(dds, geoMeans=geoMeans)
bdds=DESeq(dds, test = "Wald", fitType = "local", betaPrior=FALSE)
resultsNames(bdds)

And for contrast example.....

ctrst = colMeans(m[d$Type == “2” & d$chem1 == "Y" & d$chem2 == "N",]) -
        colMeans(m[d$Type == “2” & d$chem1 == "N" & d$chem2 == "N",])
        

res=results(bdds, contrast=ctrst)

Does this seem right? Also, I've noticed I get quite a few segfaults within R when running the deseq2/phyloseq packages... is this just unique to my setup?

[Michael Love]

There is a way to control for the animal effect in this experiment, by recoding the animals within blocks of chem1*chem2 and adding an interaction term. It gets a bit trickier to form contrasts, but I will show how to do this.

It helps to clean up the colData a bit, I'll show on a clean version, and you can reorder your samples correspondingly. Note that you should reorder the columns of the count matrix so that they match. If you have a SummarizedExperiment, the colData and matrix are tied together, so simply reordering the columns of the SummarizedExperiment does both for you.

I made a sample data.frame like your data:

d <- data.frame(anim=factor(rep(1:12,each=3)),
  chem1=factor(rep(rep(c("N","Y"),each=6),3)),
  chem2=factor(rep(rep(c("N","Y"),each=3),6)),
  tiss=factor(rep(1:3,12)))

this looks like

    anim chem1 chem2 tiss
 1     1     N     N    1
 2     1     N     N    2
 3     1     N     N    3
 4     2     N     Y    1
 5     2     N     Y    2
 6     2     N     Y    3
 7     3     Y     N    1
 8     3     Y     N    2
 9     3     Y     N    3
 10    4     Y     Y    1
 11    4     Y     Y    2
 12    4     Y     Y    3
 13    5     N     N    1
 14    5     N     N    2
 15    5     N     N    3
 16    6     N     Y    1
 17    6     N     Y    2
 18    6     N     Y    3
 19    7     Y     N    1
 20    7     Y     N    2
 21    7     Y     N    3
 22    8     Y     Y    1
 23    8     Y     Y    2
 24    8     Y     Y    3
 25    9     N     N    1
 26    9     N     N    2
 27    9     N     N    3
 28   10     N     Y    1
 29   10     N     Y    2
 30   10     N     Y    3
 31   11     Y     N    1
 32   11     Y     N    2
 33   11     Y     N    3
 34   12     Y     Y    1
 35   12     Y     Y    2
 36   12     Y     Y    3

Now reorder by chem1, chem2 and tiss. Then add a new variable which will be used to control animal specific effects within each block of chem1*chem2.

idx <- with(d, order(chem1, chem2, tiss))
d <- d[idx,]
d$anim.nested <- factor(rep(1:3,12))

You can examine d to see how this works.

The model matrix is then the chem1 and chem2 terms and interaction, as well as the interactions of these with animal, and the interaction of these with tissue:

m <- model.matrix(~ chem1 + chem2 + chem1:chem2 + chem1:chem2:anim.nested + chem1:chem2:tiss, data=d)

The column names of this are:

> colnames(m)
  [1] "(Intercept)"                "chem1Y"
  [3] "chem2Y"                     "chem1Y:chem2Y"
  [5] "chem1N:chem2N:anim.nested2" "chem1Y:chem2N:anim.nested2"
  [7] "chem1N:chem2Y:anim.nested2" "chem1Y:chem2Y:anim.nested2"
  [9] "chem1N:chem2N:anim.nested3" "chem1Y:chem2N:anim.nested3"
 [11] "chem1N:chem2Y:anim.nested3" "chem1Y:chem2Y:anim.nested3"
 [13] "chem1N:chem2N:tiss2"        "chem1Y:chem2N:tiss2"
 [15] "chem1N:chem2Y:tiss2"        "chem1Y:chem2Y:tiss2"
 [17] "chem1N:chem2N:tiss3"        "chem1Y:chem2N:tiss3"
 [19] "chem1N:chem2Y:tiss3"        "chem1Y:chem2Y:tiss3"

Many of the contrasts become hard to form, though possible with numeric contrasts. A trick for getting the numeric contrasts of the other comparisons, is to take the column means of the model matrix, for the groups of samples you want to compare. For instance, if you want to find the main chem1 effect in tissue 2, the contrast would be:

ctrst = colMeans(m[d$tiss == "2" & d$chem1 == "Y" & d$chem2 == "N",]) -
        colMeans(m[d$tiss == "2" & d$chem1 == "N" & d$chem2 == "N",])

This ends up being the chem1 base effect, the chem1Y effect in tissue 2 minus the chem1N effect in tissue 2 and a few 1/3 and -1/3 to the animal interaction effects.

Edit: note that interaction terms can be tested using the same logic, but you need to subtract the groups with only one chem and then add the base group. So it looks like: chem1Y chem2Y - chem1N chem2Y - chem1Y chem2N + chem1N chem2N.

Also, I had to update that custom line of size factor code for when all rows have a zero. If you were using that line of code to avoid the issue that all rows of your count matrix have one or more zeros, you should get the edited version:

https://support.bioconductor.org/p/62246/

We're implementing a solution here which was on the list of todos for a while now, but it won't be out for 6 months. Typical RNA-Seq doesn't have this issue, but since phyloseq we are seeing more metagenomics datasets which often have no genes without a zero.

[Michael Love]

I'll add more tomorrow, you can control for animal, but you have to reformulate the terms.

[Michael Love]

aha, I thought you meant that it was only three animals total. You can't just add a term for animal because it forms linear combinations with chem1*chem2...