[BioC] multi-level, multi-factor results interpretation.

[Ingrid Lindquist]

Hi Simon,

I am working with a multi-factor experiment that, in a sense, doesn't 
have biological replicates (replicates were done at different times, and 
the differences in the chemistry/instrumentation is enough to be a 
rather large component of variance) and has three levels.  By three 
levels, I mean there are 2 different treatments and 1 control.  When I 
push the workflow through, everything works fine, but I'm having a hard 
time delineating which level is responsible for a gene being identified 
as significantly differentially expressed. I understand that (in the 
following example) ConditionX, ConditionY, AlternativeVariableold are 
log2 fold changes in relation to "ctrl" (in the case of condition) or 
"new" (in the case of AlternativeVariable).  Can I assume that the 
largest FC of these three fields within each gene instance is likely the 
reason that gene had a padj within significance? My plan is to separate 
out the dataset so that I'm only running 2 levels at once (x vs ctrl; y 
vs ctrl), but I'm hoping you can shed light on how I can interpret these 
results I've mentioned here.

The output I'm referring to has the following headers (generated by the 
last line (write.table) in my workflow below):

Gene    Pval    Padj    Intercept     ConditionX    ConditionY    
AlternativeVariableold    Deviance Converged


My workflow:

counts <- (file, header=TRUE, row.names=1)

Design<- data.frame(
+  row.names = colnames(SampledCounts),
+ Condition = c( "x", "ctrl", "ctrl", "x", "y", "y"),
+ AlterativeVariable= c( "old", "new", "old", "new", "new", "old"))

library(DESeq)

cdsFull <- newCountDataSet (counts, Design)
cdsFull <-estimateSizeFactors(cdsFull)
cdsFull <-estimateDispersions(cdsFull, method="blind", 
sharingMode="fit-only")
fit1<-fitNbinomGLMs(cdsFull, count ~ condition + AlternativeVariable)
fit0<-fitNbinomGLMs(cdsFull, count ~ AlternativeVariable)
pvalsGLM <-nbinomGLMTest(fit1, fit0)
padjGLM <-p.adjust(pvalsGLM, method = "BH")
write.table(data.frame(geneID=row.names(counts(cdsFull)), pval=pvalsGLM, 
padj=padjGLM, fit1), file= "result.txt")

Thanks for your help,
Ingrid