[BioC] DESeq2 - paired sample, interaction and fold change

[Michael Love]

hi Samuel,


On Thu, Jun 5, 2014 at 4:35 AM, samuel collombet
<samuelcollombet at gmail.com> wrote:
>
>> yes, by including the patient variable in the design, differences in
>> counts which are due to patient differences are accounted for.
>
> I am still not sure I get it, certainly because I am not familiar with
> expanded design matrices...
> the contrast "conditionWT.treatmentCt" will define the vector of contrast
> (named c in the doc/paper) that, multiplied by the vector of GLM coefficient
> Beta, will give the vector of coefficient that will be used to calculate the
> logarithmic posterior, right?
> Is it possible to get the contrast vector / pair of contrast vector that
> correspond to a resultName?
>

For the interaction terms, like 'conditionWT.treatmentCt', there is a
single coefficient in the model, so we don't have to build a contrast.
We just pull out the coefficient and it's standard error. When we want
to contrast multiple terms, then the results() function creates a
vector 'c'. The GLM coefficient vector beta is the maximum a posterior
estimates. c * beta gives the numerator of the contrast as describes
in ?results.

If you supply a numeric contrast, the numbers should correspond to the
terms in resultsNames. This is described in ?results and there are
useful examples there.


>>
>>> Is it possible to get the individual fold change for each patients?
>>
>> Can you be more specific on which fold change? Are you interested in the
>> tumor-treatment interaction term for each patient?
>
> yes
>
>>   How many replicates do you have for each combination of patient x
>> condition x treatment?
>
> 3
>
>>   If not many, I'd recommend sticking to the general effect above.
>
> We suspect that not all sample respond to the treatment, and so we would
> like to know if the 'non-significance' of some genes is due to global
> non-response, or 1 sample that is not responding at all (or if they are all
> responding differently, so the variability is too high to be significant).
> I can get this kind of information by getting the rlogTransformation of the
> counts and comparing the different samples, but these will not be the exact
> values used for the test... I would like to be able to extract these last
> ones.
>

If you have replicates for each combination of patient, treatment and
condition (I'd recommend at least 3), then you can use a model with
2nd order interactions. You will have to turn off the betaPrior, as we
do not now support shrinkage of log fold changes for models with 2nd
or higher order interactions.

Here is some example code:

# construct an example dataset with 8 combinations x 3 replicates each
dds <- makeExampleDESeqDataSet(n = 1000, m = 24)
dds$patient <- factor(rep(c("1","2"),each=12),levels=c("1","2"))
dds$condition <-
factor(rep(rep(c("WT","tumor"),each=6),2),levels=c("WT","tumor"))
dds$treatment <- factor(rep(rep(c("Ct","A"),each=3),4),levels=c("Ct","A"))
as.data.frame(colData(dds))

# allow interactions between all factors
design(dds) <- ~ patient*condition*treatment

# turn off beta prior
dds <- DESeq(dds, betaPrior=FALSE)

# just to visualize the model matrix in this case
(mm <- unname(attr(dds,"modelMatrix")))
par(mar=c(15,1,1,1))
image(t(mm[nrow(mm):1,]),xaxt="n",yaxt="n")
axis(1,0:7/7,resultsNames(dds),las=2)

# the results names
resultsNames(dds)

# this is the interaction term for patient 1
# i.e. does treatment A have a different effect in tumor for patient 1
results(dds, name="conditiontumor.treatmentA")

# this is the interaction term for patient 2
# i.e. does treatment A have a different effect in tumor for patient 2
# note: we add the interaction from before with the additional
interaction for patient 2
# see ?results for details on providing the list()
results(dds, contrast=list(c("conditiontumor.treatmentA",
 "patient2.conditiontumor.treatmentA"),character(0)))

# the previous contrast is the same as constructing the following
numeric contrast
ctrst <- as.numeric(resultsNames(dds) %in% c("conditiontumor.treatmentA",
 "patient2.conditiontumor.treatmentA"))
ctrst
results(dds, contrast=ctrst)


Mike


> Many thanks!