Hi all,

We are running into a problem using an ANODEV with DESeq.

Let me state our data, our aim, and our dilemma.

Our data (time course with various times titled A, B and C):

A design matrix would look something like this (this is a dramatic simplification, our real data contrains over 200 participants) :

Design <- data.frame(row.names=colnames(CountTable)
, timecourse = c("A","A","A","A","A","A",
"B","B","B","B","B","B","C","C","C","C","C","C")
, pairs=c("1", "2", "3", "4", "5", "6","1", "2", "3", "4",
"5","6","1", "2",
"3", "4", "5", "6")
, condition=c("Control", "Control", "Control", "Disease", "Disease",
"Disease","Control", "Control", "Control", "Disease", "Disease",
"Disease","Control", "Control", "Control", "Disease", "Disease",
"Disease"))
countTable <- CountTable
cds <- newCountDataSet(countTable, Design)
cds <- estimateSizeFactors(cds)
sizeFactors( cds )
cds <- estimateDispersions(cds, method="pooled-CR", fitType="local")

Our aim:

What we aim to do is to test for DE of transcripts across all 3 time points for disease and controls seperatly (using DESeq ANODEV) but we want to be able to identify at which time points these transcripts are being DE. In other words, we want to compare DE transcripts with respect to specific time points between cases and controls. Our remaining code looks like this:

fit0 <- fitNbinomGLMs (cds, count ~ timecourse)
fit1 <- fitNbinomGLMs ( cds, count ~ timecourse + condition )
str(fit1)
pvals <- nbinomGLMTest( fit1, fit0 )
padj <- p.adjust( pvals, method="fdr" )
rownames(counts(cds))[ padj < .1 ]

Our delimma:

ANODEV is in fact yielding a list of significantly differentially expressed genes. However, it is not informative if we want to know from which time point this differential expression occurred. From our knowledge there is no implementation of post-hoc testing(or equivilent) after ANODEV in DESeq that states at which time point a gene is being called significant from within a timecourse study.

What insight can anyone offer us regarding DESeq towards better understanding at which time points are our DEGs are related to? In short, how can one tell at which time points (or both) is DESeq calling a transcript significantly DE?

Also, does this approach take into consideration the matching paired data?

Yours,
Michael

Dear Michael,

What you want to do is easy and fast using the edgeR package, without any need for ad hoc workarounds like subsetting your data. See McCarthy et all (NAR 2012):

  http://www.ncbi.nlm.nih.gov/pubmed/22287627

Your experimental design seems to fit exactly into the framework of Section 3.3 of the edgeR User's Guide.

Best wishes

Gordon

Gordon!

On second pass, edgeR doesn't seem to address my question.

We have disease and control data spread over 3 time points (baseline-T0, time after disease causing event-T1, elapsed time after event-T3) for each of these two conditions. I want to be able to find DEGs accordingly to condition across time-points to decipher which times (either at treatment or elapsed time after treatment) our genes are DE.

Any further insight would be fantastic.

Yours,
Michael