Dear all, We are currently analyzing RNA-seq data using DEseq and would like to ask for your advice regarding the possibility to treat a variable as a "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)* *Diet_effect_adjusted< p.adjust(Diet_effect, method="BH" )* 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)* *Interaction_effect _adjusted <- p.adjust(Interaction_effect, method="BH" )* 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 -- Akiko Sugio, Ph.D. INRA, UMR Institut de Génétique, Environnement et Protection des Plantes (IGEPP) Domaine de la Motte, 35653 Le Rheu cedex - France tel: 00 33 (0) 2 23 48 51 53 akiko.sugio@rennes.inra.fr [[alternative HTML version deleted]]

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 > ask for your advice regarding the possibility to treat a variable as a > "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)* > > *Diet_effect_adjusted< p.adjust(Diet_effect, method="BH" )* > > 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)* > > *Interaction_effect _adjusted <- p.adjust(Interaction_effect, method="BH" )* > > 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 >