[BioC] how to design matrix by edgeR

[Gordon K Smyth]

Dear Shao Gao,

I don't really understand your questions, and I'm afraid I have no more 
time to give you tutorials at the moment.

I will just say that the formula ~time+treatment (that you did not mention 
in earlier emails) estimates the treatment effect averaged over all the 
different times.  It estimates just one average treatment effect.

The formula I suggested, ~time+time:treatment estimates a separate 
treatment effect for each time.  In other words, it estimates seven 
different time-specific treatment effects, one for 0h, one for 6h, etc. 
It seems clear from the experimental design that this is what is required.

If you find that using model formula in R is too obscure to understand, 
then here is a longer method that you might find more intuitive.  It 
treats your experiment as a oneway layout, from which you can extract all 
the comparison you want.  You create a combined factor to label all your 
samples.  Lets say you use "T16" to mean treatment at 16h and C16 to mean 
control at 16h, and so on.  So create a combined factor with something 
like:

   Group <- factor(c("T0","T0",...,"C48"))

Then

   Group <- relevel(Group, ref="T0")
   design <- model.matrix(~0+Group)
   colnames(design) <- levels(Group)

Fit the full model (after estimating the dispersions):

   fit <- glmFit(y,design)

Then later on you can make any comparison you want.  Suppose you want test 
for a treatment effect at time 16, that is T16-C16.  You can ask for this 
by:

   T16vsC16 <- makeContrasts(T16-C16, levels=design)
   lrt <- glmFit(y,fit,contrast=as.vector(T16vsC16))

This gives exactly the same results as the approach I suggested in my last 
email.

You can use this for any comparison at all.  For example, is the treatment 
effect larger at 6h than at 0h?

   con <- makeContrasts((T6-C6)-(T0-C0), levels=design)
   lrt <- glmFit(y,fit,contrast=as.vector(con))

I hope this gives you enough to go on with.  If not, I hope that others 
can help you.

Best wishes
Gordon

---------------------------------------------
Professor Gordon K Smyth,
Bioinformatics Division,
Walter and Eliza Hall Institute of Medical Research,
1G Royal Parade, Parkville, Vic 3052, Australia.
Tel: (03) 9345 2326, Fax (03) 9347 0852,
smyth at wehi.edu.au
http://www.wehi.edu.au
http://www.statsci.org/smyth

On Tue, 22 May 2012, wang peter wrote:

> dear Dr.Gordon:
>   thank  you for your patient reply, but i still donot understand
> what problem in my factor design?i have 3 questions as below:
>
> the 1st question:
>   my experiments involves two factors, time and treatment.
>   so i define:
> treatment=factor(c(rep('treated',24),rep('control',11)))
> time=factor( c('0h','0h','0h','6h','6h','6h','6h','12h',
> '12h','12h','12h','18h','18h','18h','18h','24h','24h','24h','36h','36h','36h','48h','48h','48h','0h','0h','0h','6h','12h','18h','24h','36h','48h'))
>
> the 2nd question:
> i want to get DE cross time considersing control.
> and my model.matrix is
> design <- model.matrix(~time+treatment)?
> or design <- model.matrix(~treatment+time)?
>
> and yours is design <- model.matrix(~time+time:treatment)
> which one is correct?
>
>
> the third question:
>      i donot do what is the difference between test by glmLRT
> and test by pairedwised test functions?
>
>
> To test for treatment effect at time 0:
>
>  lrt0 <- glmLRT(y,fit,coef=8)
>  topTags(lrt0)
>
> To test for treatment effect at time 1:
>
>  lrt1 <- glmLRT(y,fit,coef=9)
>  topTags(lrt1)
>
> thank you very much
> if possible, can you give me some detailed material
> of glmLRT or glmFIT
> now i can only read your paper published recently
>
> shan gao
>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}