[BioC] mixed-effects analysis within DEseq

Mike Love love at molgen.mpg.de
Mon Oct 8 20:19:14 CEST 2012

Hi Akiko,

I believe that DESeq treats all effects as "fixed", i.e. parameters to
be estimated.

Though with such small number of groups and replicates, it should not
make a big difference in the estimates or the testing of the treatment
effect.  For example, below I fit a generalized linear model and a
generalized linear mixed model to Poisson counts of the design that you
have (2 x 3 x 2).  The coefficient for the fixed effect, x, is almost
identical, as is the change in deviance between the full and reduced
model (this 435.33 which is the basis for the Chi-squared test).  I did
not explore what might happen with interaction terms.

> n <- 12
> x <- factor(rep(0:1,each=n/2))
> z <- factor(rep(c(0:2),times=n/3))
> y <- rpois(n, exp(1 + as.integer(x) + as.integer(z)))
> glm.fit1 <- glm(y ~ x + z, family=poisson)
> coef(glm.fit1)
(Intercept)          x1          z1          z2
2.919019    1.105920    0.949297    2.009069

> library(lme4)
> glmix.fit1 <- glmer(y ~ x + (1|z), family=poisson)
> coef(glmix.fit1)
\$z
(Intercept)       x1
0    2.928748 1.105916
1    3.868466 1.105916
2    4.926723 1.105916

> glm.fit0 <- glm(y ~ z, family=poisson)
> glmix.fit0 <- glmer(y ~ (1|z), family=poisson)
> anova(glm.fit0, glm.fit1, test="Chi")
Analysis of Deviance Table

Model 1: y ~ z
Model 2: y ~ x + z
Resid. Df Resid. Dev Df Deviance  Pr(>Chi)
1         9     448.91
2         8      13.58  1   435.33 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> anova(glmix.fit0, glmix.fit1)
Data:
Models:
glmix.fit0: y ~ (1 | z)
glmix.fit1: y ~ x + (1 | z)
Df    AIC    BIC   logLik  Chisq Chi Df Pr(>Chisq)
glmix.fit0  2 472.70 473.67 -234.351
glmix.fit1  3  39.37  40.82  -16.685 435.33      1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

best,

Mike

On 10/6/12 12:00 PM, bioconductor-request at r-project.org wrote:
>
> Message: 4
> Date: Fri, 05 Oct 2012 16:15:18 +0200
> From: Akiko Sugio <akiko.sugio at rennes.inra.fr>
> To: bioconductor at r-project.org
> Subject: [BioC] mixed-effects analysis within DEseq
> Message-ID: <506EEB76.1020102 at rennes.inra.fr>
> Content-Type: text/plain
>
> Dear all,
>
> We are currently analyzing RNA-seq data using DEseq and would like to
> "random effect" rather than a "fixed effect".
>
> In our experiment, we are working with clonal individuals (aphids). We
> are working with three different clones. Each of these clones is fed
> with two different diets (treatment).We have two biological replicates
> for each condition. Hence, we have a total of 12 RNAseq samples (3
> different clones x two diet treatment x 2 replicates). We are mainly
> interested in identifying genes that are differentially expressed in
> response to the diet.
>
> Currently, we are running GLM models in DESeq:
>
> *Model0 <- fitNbinomGLMs( cds, count ~diet * cloneID )*
>
> We then test for the significance of each variable (i.e. cloneID, diet,
> and the interaction diet:cloneID) by comparing the model including the
> variable of interest to the model estimated without this variable.
>
> To test for the effect of the diet:
>
> *Model1 <- fitNbinomGLMs( cds, count ~cloneID)*
>
> *Model2<-fitNbinomGLMs( cds, count ~cloneID + diet )*
>
> *Diet_effect <- nbinomGLMTest(Model2, Model1)*
>
>
> To test for the effect of the interaction:
>
> *Model0 <- fitNbinomGLMs( cds, count ~diet * cloneID ).*
>
> *Model3 <- fitNbinomGLMs( cds, count ~diet + cloneID ).*
>
> *Interaction_effect<- nbinomGLMTest(Model0, Model3)*
>
>
> For genes with a significant interaction between diet and CloneID, we do
> not interpret the effect of each variable (i.e. diet and CloneID).
>
> We are however not entirely satisfied by our model, because we think it
> would be better to treat the variable cloneID as a random effect. Is it
> possible to conduct such mixed-effects analysis within DEseq?
>
> Thank you very much in advance.
> Akiko
>