[BioC] edge R design matrix pairwise comparison at different time points after infection with replicates
Gordon K Smyth
smyth at wehi.EDU.AU
Wed Jun 27 00:57:58 CEST 2012
On Tue, 26 Jun 2012, Kaat De Cremer wrote:
> Dear Gordon,
>
> I do normalize the data, after filtering for low expressed genes. Below
> you can see my first y$samples (when I forgot to correct the library
> size) and the corrected version. This correction for library size (which
> is small and of the same order for each library) causes a big difference
> in number of differentially expressed genes. But what I find odd is that
> I don't see this difference when I analyze each time point separate,
> including only data of one time point in edgeR. I notice that the
> normalization factors are more consistent within each time point
> (especially for the first time point), could this cause the difference
> between the different approaches?
No.
> group lib.size norm.factors
> 12hpi C1 C 9498975 1.2036962
> 12hpi C2 C 5757034 1.1233049
> 12hpi C3 C 6367418 1.1285791
> 12hpi T1 T 6250384 1.1524029
> 12hpi T2 T 5606716 1.1351563
> 12hpi T3 T 7294081 1.1001216
> 24hpi C1 C 5873429 0.9832233
> 24hpi C2 C 6357801 0.8988242
> 24hpi C3 C 6778635 0.9467255
> 24hpi T1 T 6596149 1.0076682
> 24hpi T2 T 5451307 0.8962813
> 24hpi T3 T 5726377 0.9136752
> 48hpi C1 C 8232847 0.9945134
> 48hpi C2 C 8079906 0.9643586
> 48hpi C3 C 7632570 0.9228554
> 48hpi T1 T 15413395 0.9086659
> 48hpi T2 T 6289516 0.9171062
> 48hpi T3 T 5003424 0.8942136
>
> Corrected for library size:
>
> group lib.size norm.factors
> 12hpi C1 C 9485607 1.2036962
> 12hpi C2 C 5749033 1.1233049
> 12hpi C3 C 6358791 1.1285791
> 12hpi T1 T 6242328 1.1524029
> 12hpi T2 T 5598562 1.1351563
> 12hpi T3 T 7285598 1.1001216
> 24hpi C1 C 5866969 0.9832233
> 24hpi C2 C 6351387 0.8988242
> 24hpi C3 C 6771865 0.9467255
> 24hpi T1 T 6588283 1.0076682
> 24hpi T2 T 5445788 0.8962813
> 24hpi T3 T 5721236 0.9136752
> 48hpi C1 C 8223672 0.9945134
> 48hpi C2 C 8071573 0.9643586
> 48hpi C3 C 7624662 0.9228554
> 48hpi T1 T 15393879 0.9086659
> 48hpi T2 T 6282481 0.9171062
> 48hpi T3 T 4997554 0.8942136
>
>
> About the different approaches of analyzing the data (all data together
> or each time point separate), I included my MDS plot in attachment. I do
> see more variability at the later time point (especially between treated
> and control samples).
This isn't variability. It's a treatment effect. The treated samples at
48h are closely clustered showing low variability.
The MDS plots hows that you have a treatment effect mainly at time 48h.
> To estimate the dispersion, does edgeR compare all different samples as
> if they were independent, or does edgeR bring in account which samples
> are of the same time point and same treatment?
Yes, certainly it takes into account which samples are of the same time
and the same treatment. It uses the design matrix you give it.
> still feel that my 3 time points look very different and that combining
> all data will overestimate the dispersion.
Of course the time points are different, but because of systematic effects
not because of random variability.
Best wishes
Gordon
> Thank you so much for your answers,
>
> Kaat
>
>
>
> -----Original Message-----
> From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU]
> Sent: dinsdag 26 juni 2012 2:04
> To: Kaat De Cremer
> Cc: Mark Robinson; Bioconductor mailing list
> Subject: RE: design matrix edge R pairwise comparison at different time points after infection with replicates
>
> Dear Kaat,
>
> 1) It is not generally true that you will find more DE genes analysing
> just one time separately rather than using all the libraries in one
> linear model. It is possible in principle that the later time may show
> more variability, and that might justify separate analysis of the
> different times. However I would do that only when you have a good a
> priori biological reason for expecting such an effect and data
> exploration (such as an MDS plot) confirms it. Otherwise the extra
> instability of dispersion estimation with a small number of libraries is
> not justified. I would not make such decisions merely on the basis of
> which analysis gives more DE genes.
>
> 2) I don't think there is yet a definitive rule regarding filtering and
> library sizes, but I prefer to recompute library sizes and scale
> normalize (calcNormFactors) after filtering. Recomputing the library
> sizes doesn't make a lot of difference, because scale normalization will
> self correct anyway.
>
> Are you normalizing your data?
>
> Best wishes
> Gordon
>
> On Mon, 25 Jun 2012, Kaat De Cremer wrote:
>
>> Again,
>> Thank you both very much for your reply.
>>
>> I have analyzed my time course data now in several different ways and
>> noticed some differences:
>>
>>
>> 1) when I analyze the data of all time points together and look for DE
>> genes at one time point, I find less DE genes compared to when I use
>> only the data of that one time point in edgeR. I assume this is
>> because the dispersion is larger when I include all the different time
>> points at once? In that case, is this the right way to go? The
>> dispersion can be larger at one time point compared to another I would think.
>>
>> 2) I noticed in the edgeR user's guide that in some examples you
>> correct the library size after filtering for low expressed genes, and
>> in other examples you don't. Correcting this library size gives less
>> DE genes for my data at all time points when I analyze all data
>> together and then look for DE genes at each time point. I don't see
>> this difference when I only look at the data of one time point in edgeR.
>>
>> I hope you can comment on this,
>>
>>
>> Thank you,
>> Kaat
>>
>>
>>
>> -----Original Message-----
>> From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU]
>> Sent: zondag 24 juni 2012 2:42
>> To: Kaat De Cremer
>> Cc: Bioconductor mailing list; Mark Robinson
>> Subject: design matrix edge R pairwise comparison at different time
>> points after infection with replicates
>>
>> Hi Kaat,
>>
>> I'll jump in and continue on from Mark's help.
>>
>> To test for treatment effects separately at each time, the easiest way is to include the terms "time+time:treat" in your model formula.
>>
>> I'll assume that your "tests" are independent replicates of the whole experiment. If there are batch effects associated with the tests that you need to correct for, then your complete design matrix might be:
>>
>> design <- model.matrix(~test+time+time:treat)
>>
>> This produces a design matrix with the following columns:
>>
>> > colnames(design)
>> [1] "(Intercept)" "test2" "test3" "time24hpi"
>> [5] "time48hpi" "time12hpi:treatT" "time24hpi:treatT" "time48hpi:treatT"
>>
>> So testing for treatment effects at each time is easy. To test for treatment effect as time 12h:
>>
>> fit <- glmFit(y, design)
>> lrt <- glmLRT(y, fit, coef="time12hpi:treatT")
>>
>> etc. To test for treatment effect at time 24h:
>>
>> lrt <- glmLRT(y, fit, coef="time24hpi:treatT")
>>
>> and so on.
>>
>> Best wishes
>> Gordon
>>
>>> Date: Fri, 22 Jun 2012 13:11:41 +0000
>>> From: Kaat De Cremer <Kaat.DeCremer at biw.kuleuven.be>
>>> To: Mark Robinson <mark.robinson at imls.uzh.ch>
>>> Cc: bioconductor list <bioconductor at r-project.org>
>>> Subject: Re: [BioC] design matrix edge R pairwise comparison at
>>> different time points after infection with replicates
>>>
>>> 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","48h
>>> p
>>> i","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
>>>>
>>
>> ______________________________________________________________________
>> The information in this email is confidential and intended solely for the addressee.
>> You must not disclose, forward, print or use it without the permission of the sender.
>> ______________________________________________________________________
>>
>
> ______________________________________________________________________
> The information in this email is confidential and inte...{{dropped:10}}
More information about the Bioconductor
mailing list