[BioC] design matrix edge R pairwise comparison at different time points after infection with replicates
Kaat De Cremer
Kaat.DeCremer at biw.kuleuven.be
Fri Jun 22 15:11:41 CEST 2012
Hi Mark,
Thank you for your suggestion,
I really appreciate your time.
Working in R is new to me so it has been a struggle using edgeR, but I think I managed it using only 2 factors (test and treatment). Now that I will be including 3 factors (test, treatment and time) in one analysis it is clear to me that I still don't understand how it works exactly.
Below you can see my workspace with the only design matrix I could come up with, but I don't see which coefficients I should include or which contrast vector to use in the glmLRT function to make the comparison of control-treatment at each time point separate, ignoring the other 2 time points. Is this possible with this design matrix? Or is the matrix wrong for this purpose?
Thanks!
Kaat
> head(x)
12hpi C1 12hpi C2 12hpi C3 12hpi T1 12hpi T2 12hpi T3 24hpi C1 24hpi C2
Lsa000001.1 0 1 1 2 0 2 1 1
Lsa000002.1 5 4 0 5 6 6 6 4
Lsa000003.1 10 9 7 5 5 8 6 2
Lsa000004.1 1 1 1 1 1 1 1 3
Lsa000005.1 1 0 1 0 2 0 0 1
Lsa000006.1 510 223 228 287 222 268 303 358
24hpi C3 24hpi T1 24hpi T2 24hpi T3 48hpi C1 48hpi C2 48hpi C3 48hpi T1
Lsa000001.1 0 1 1 0 0 0 0 2
Lsa000002.1 7 5 2 5 10 6 12 12
Lsa000003.1 7 5 4 2 6 5 8 2
Lsa000004.1 1 3 1 2 1 3 2 3
Lsa000005.1 0 1 0 0 1 0 0 2
Lsa000006.1 372 362 237 320 472 440 411 858
48hpi T2 48hpi T3
Lsa000001.1 0 0
Lsa000002.1 1 5
Lsa000003.1 1 0
Lsa000004.1 0 2
Lsa000005.1 1 0
Lsa000006.1 375 275
> treat<-factor(c("C","C","C","T","T","T","C","C","C","T","T","T","C","C","C","T","T","T"))
> test<-factor(c(1,1,2,3,1,2,3,2,3,1,2,3,1,2,3,1,2,3))
time<-factor(c("12hpi","12hpi","12hpi","12hpi","12hpi","12hpi","24hpi","24hpi","24hpi","24hpi","24hpi","24hpi","48hpi","48hpi","48hpi","48hpi","48hpi","48hpi"))
> y<-DGEList(counts=x,group=treat)
Calculating library sizes from column totals.
> cpm.y<-cpm(y)
> y<-y[rowSums(cpm.y>2)>=3,]
> y<-calcNormFactors(y)
design<-model.matrix(~test+treat+time)
> design
(Intercept) test2 test3 treatT time24hpi time48hpi
1 1 0 0 0 0 0
2 1 1 0 0 0 0
3 1 0 1 0 0 0
4 1 0 0 1 0 0
5 1 1 0 1 0 0
6 1 0 1 1 0 0
7 1 0 0 0 1 0
8 1 1 0 0 1 0
9 1 0 1 0 1 0
10 1 0 0 1 1 0
11 1 1 0 1 1 0
12 1 0 1 1 1 0
13 1 0 0 0 0 1
14 1 1 0 0 0 1
15 1 0 1 0 0 1
16 1 0 0 1 0 1
17 1 1 0 1 0 1
18 1 0 1 1 0 1
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$test
[1] "contr.treatment"
attr(,"contrasts")$treat
[1] "contr.treatment"
attr(,"contrasts")$time
[1] "contr.treatment"
> y<-estimateGLMCommonDisp(y,design,verbose=TRUE)
Disp = 0.07299 , BCV = 0.2702
> y<-estimateGLMTrendedDisp(y,design)
Loading required package: splines
> y<-estimateGLMTagwiseDisp(y,design)
Warning message:
In maximizeInterpolant(spline.pts, apl.smooth[j, ]) :
max iterations exceeded
> fit<-glmFit(y,design)
-----Original Message-----
From: Mark Robinson [mailto:mark.robinson at imls.uzh.ch]
Sent: vrijdag 22 juni 2012 12:03
To: Kaat De Cremer
Cc: bioconductor list
Subject: Re: [BioC] design matrix edge R pairwise comparison at different time points after infection with replicates
Hi Kaat,
It is probably better to fit all your data with a single call to glmFit(), over all 18 samples; you can test the differences of interest trough the 'coef' or 'contrast' argument on glmLRT(). That would afford you more degrees of freedom and presumably better estimates of dispersion, and so on.
>From your description, I can't quite figure out your design matrix. You have three factors: treatment, test and time point. First, you need to input all 18 samples and extend your 'treatment' and 'test' factor variables to have 18 values (corresponding to the columns of your table). And, then also include a time variable in your design. Some decisions might need to be made about interactions to include.
Hope that gets you started.
Best,
Mark
----------
Prof. Dr. Mark Robinson
Bioinformatics
Institute of Molecular Life Sciences
University of Zurich
Winterthurerstrasse 190
8057 Zurich
Switzerland
v: +41 44 635 4848
f: +41 44 635 6898
e: mark.robinson at imls.uzh.ch
o: Y11-J-16
w: http://tiny.cc/mrobin
----------
http://www.fgcz.ch/Bioconductor2012
On 21.06.2012, at 11:42, Kaat De Cremer wrote:
> Dear all,
>
>
> I am using edgeR to find genes differentially expressed between infected and mock-infected control plants, at 3 time points after infection.
> I have RNAseq data for 3 independent tests, so for every single test I have 6 libraries (control + infected at 3 time points).
> Having three replicates this makes 18 libraries in total.
>
> What I did until now is look at each time point separate and calculate DEgenes at that time point as shown in this script:
>
>> head(x)
> C1 C2 C3 T1 T2 T3
> 1 0 1 2 0 0 0
> 2 13 6 4 10 8 12
> 3 17 16 9 10 8 11
> 4 2 1 2 2 3 2
> 5. 1 3 1 2 1 3 0
> 6 958 457 438 565 429 518
>
>> treatment<-factor(c("C","C","C","T","T","T"))
>> test<-factor(c(1,2,3,1,2,3))
>> y<-DGEList(counts=x,group=treatment)
> Calculating library sizes from column totals.
>> cpm.y<-cpm(y)
>> y<-y[rowSums(cpm.y>2)>=3,]
>> y<-calcNormFactors(y)
>> design<-model.matrix(~test+treat)
>> y<-estimateGLMCommonDisp(y,design,verbose=TRUE)
> Disp = 0.0265 , BCV = 0.1628
>> y<-estimateGLMTrendedDisp(y,design)
> Loading required package: splines
>> y<-estimateGLMTagwiseDisp(y,design)
>> fit<-glmFit(y,design)
>> lrt<-glmLRT(y,fit)
>
>
> This works fine but I wonder if I should do the analysis of the different time points all at once? Will this make a difference?
> Unfortunately I cannot figure out how to design the matrix.
>
> I hope someone can help me,
>
> Kaat
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list