DESeq2:could not find function "Error"

[learnholickyle]

I am using DESeq2 to analysis a longitudinal data. we have two species for two treatment (control and treatment). we collect data on 0, 1, 2, 5 days three replicates. so we have 2 * 2 * 3=12 subjects' data on four time points.

Using anova, we design the formula, like this: count~species*treatment*time+Error(subject)
It works for aov. However, DESeq2 returns: Error in eval(expr, envir, enclos) : could not find function "Error"

Anyone have an idea how to solve this problem?

Thanks for your time!

kyle

[dpryan]

DESeq2 doesn't do mixed-effect models. The closest you can get would be:

count~species*treatment*time+subject

You might also consider using limma, specifically the duplicateCorrelation() function, which is often used in these sorts of cases.

[learnholickyle]

DESeq2 doesn't do mixed-effect models. The closest you can get would be:

count~species*treatment*time+subject

You might also consider using limma, specifically the duplicateCorrelation() function, which is often used in these sorts of cases.

Thanks for your help.

I am trying to use count~species*treatment*time+treatment:subject. I will try to use limma later.

[Michael Love]

hi,

DESeq2 offers fixed effects modeling, where a single dispersion value is estimated per gene. This gives the expected variance of the counts around the expected value for each sample as Var(K_ij) = mu_ij + alpha_i * mu_ij^2, where K_ij is the count for the i-th gene and the j-th sample, and alpha_i is the final estimate of dispersion produced by the estimateDispersions() function. So the variance depends on the expected value, mu_ij, but we do not estimate different dispersions for different samples.

What is your biological question of interest here?

Sometimes there are more powerful approaches than looking at all interactions between all variables, e.g., species*treatment*time.

Mike

[learnholickyle]

hi,

DESeq2 offers fixed effects modeling, where a single dispersion value is estimated per gene. This gives the expected variance of the counts around the expected value for each sample as Var(K_ij) = mu_ij + alpha_i * mu_ij^2, where K_ij is the count for the i-th gene and the j-th sample, and alpha_i is the final estimate of dispersion produced by the estimateDispersions() function. So the variance depends on the expected value, mu_ij, but we do not estimate different dispersions for different samples.

What is your biological question of interest here?

Sometimes there are more powerful approaches than looking at all interactions between all variables, e.g., species*treatment*time.

Mike

Hi, Mike:

Thanks for your reminding. we want to find the DE genes of the tolerant species (one is tolerant, the other is not) under treatment. we tried different time points, want to find more these genes.

kyle

[Michael Love]

at first glance, I would start with designs of the form

~ subject + time + species + condition+ species:condition

where the base level of species is the intolerant species, and the base level of condition (i'm calling this variable "condition" for clarity) is control.

Then run:

dds <- DESeq(dds, test="LRT", reduced = ~ subject + time + species + condition)

which will perform a likelihood ratio test of the significance of the interaction: (species = tolerant & condition = treatment).

This will produce low p-values for genes in which the tolerant species has different log fold change of expression of treatment vs control compared to the intolerant species.

[learnholickyle]

at first glance, I would start with designs of the form

~ subject + time + species + condition+ species:condition

where the base level of species is the intolerant species, and the base level of condition (i'm calling this variable "condition" for clarity) is control.

Then run:

dds <- DESeq(dds, test="LRT", reduced = ~ subject + time + species + condition)

which will perform a likelihood ratio test of the significance of the interaction: (species = tolerant & condition = treatment).

This will produce low p-values for genes in which the tolerant species has different log fold change of expression of treatment vs control compared to the intolerant species.

Thanks Mike! That is great helpful.

I am reading the manual about DESeq2 you wrote on Feb. 19 2014. Could I ask one question?

According to the PCA analysis (on Page 18), "condition" seems to contribute a lot to PC1. While after the likelihood ratio test (on Page 22), there is no big change between 'type+condition' and 'type' (only 4 more genes after LRT). It means 'condition' is not important.
There must be sth I did not understand in these two places.

Thanks for your help!

kyle

[dpryan]

In addition to Michael's most recent reply, you might also consider excluding timepoint 0. Unless there's an acute reaction to the treatment in the sensitive group, this is unlikely to be informative.

Regarding the results that you're seeing on page 22 of the vignette (I'm obviously not Michael, but I'll reply anyway), that doesn't mean that 'condition' is unimportant, just that the results from using a Wald test and an LRT (at least in terms of p-values being <0.1 or not) are relatively similar. There are >1000 DE genes due to condition in both cases, so I'd argue that 'condition' is probably quite important.

[learnholickyle]

In addition to Michael's most recent reply, you might also consider excluding timepoint 0. Unless there's an acute reaction to the treatment in the sensitive group, this is unlikely to be informative.

Regarding the results that you're seeing on page 22 of the vignette (I'm obviously not Michael, but I'll reply anyway), that doesn't mean that 'condition' is unimportant, just that the results from using a Wald test and an LRT (at least in terms of p-values being <0.1 or not) are relatively similar. There are >1000 DE genes due to condition in both cases, so I'd argue that 'condition' is probably quite important.

Hi Devon:

Thanks for your explanation of LRT results. It is similar, not unimportant. Good points!

kyle

[Michael Love]

Yes, Devon's got it right.