[BioC] how to design matrix by edgeR
Gordon K Smyth
smyth at wehi.EDU.AU
Wed May 23 01:41:54 CEST 2012
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}}
More information about the Bioconductor
mailing list