[BioC] paired samples and time course experiments in LIMMA

James W. MacDonald jmacdon at uw.edu
Fri Oct 25 15:56:34 CEST 2013


Hi Lauren,

On Thursday, October 24, 2013 11:53:12 PM, Lauren Sassoubre wrote:
> Hello,
> I have a question about which analysis to do in LIMMA. I did three experiments (three biological replicates) each with one treatment microcosm ("L") and one control microcosm ("D"). For each experiment, I sampled both the treatment and control microcosms at 5 time points (0 hr, 2 hr, 6 hr, 12 hr, 24 hr). With 3 replicate experiments * (1 treatment + 1 control) * 5 time points, I have data from 30 microarrays.
>
> Based on section 8.6 of the LIMMA manual (which is incredibly helpful thank you!) here is the way I set up my targets table:
> FileName	Target
> EF731D0.CEL	D.0hr
> EF731D2.CEL	D.2hr
> EF731D6.CEL	D.6hr
> EF731D12.CEL	D.12hr
> EF731D24.CEL	D.24hr
> EF731L0.CEL	L.0hr
> EF731L2.CEL	L.2hr
> EF731L6.CEL	L.6hr
> EF731L12.CEL	L.12hr
> EF731L24.CEL	L.24hr
> EF813D0.CEL	D.0hr
> EF813D2.CEL	D.2hr
> EF813D6.CEL	D.6hr
> EF813D12.CEL	D.12hr
> EF813D24.CEL	D.24hr
> EF813L0.CEL	L.0hr
> EF813L2.CEL	L.2hr
> EF813L6.CEL	L.6hr
> EF813L12.CEL	L.12hr
> EF813L24.CEL	L.24hr
> EF815D0.CEL	D.0hr
> EF815D2.CEL	D.2hr
> EF815D6.CEL	D.6hr
> EF815D12.CEL	D.12hr
> EF815D24.CEL	D.24hr
> EF815L0.CEL	L.0hr
> EF815L2.CEL	L.2hr
> EF815L6.CEL	L.6hr
> EF815L12.CEL	L.12hr
> EF815L24.CEL	L.24hr
>
> I would like to determine which genes changed over time in the treatment relative to the control at each time point.  Is the following code correct to answer this question??

Yes, if I understand your question correctly. These are all interaction 
terms, where you are looking for genes that respond differently over 
time to treatment as compared to control.

>
> cont.dif=makeContrasts(Dif2hr=(L.2hr-L.0hr)-(D.2hr-D.0hr), Dif6hr=(L.6hr-L.0hr)-(D.6hr-D.0hr), Dif12hr=(L.12hr-L.0hr)-(D.12hr-D.0hr), levels=design)
> fitdif=contrasts.fit(fit, cont.dif)
> fitdif=eBayes(fitdif)
> topTable(fitdif, genelist=genetext, adjust="BH")
>
> I also would like to determine which genes were differentially expressed between the treatment (L) and control (D) at each time point but I'm not sure which of the following two options to use??
>
> option #1:
> contLD.matrix=makeContrasts(L.0hr-D.0hr, L.2hr-D.2hr, L.6hr-D.6hr, L.12hr-D.12hr, L.24hr-D.24hr, levels=design)
> fitLDcontmatrix=contrasts.fit(fit, contLD.matrix)
> fitLDcontmatrixB=eBayes(fitLDcontmatrix)
>
> OR
>
> option #2: Since each biological replicate experiment starts with the same bacterial culture (half exposed to the treatment while half are kept as a control), I started to think that the treatment and control samples from each experiment are paired. This made me think about doing paired sample tests for each time point. So, following the paired sample section of the LIMMA manual I created a different targets table and did the following code for each time point:
>

With cell culture, you are somewhere in between biological and 
technical replicates, regardless. In other words each cell culture is 
in its own milieu, so hypothetically you are introducing some sort of 
biological variability, but that is far different from one would 
usually consider a true biological replicate.

As with all statistical analyses the first question you have to answer 
is what assumptions you are willing to make, and why? You could blindly 
assume that the samples are highly correlated, and pair them, or you 
could assume that they are not that correlated, and assume 
independence. Or you could use the duplicateCorrelation function in 
limma to estimate the correlation structure, and determine from that 
what assumptions you think are reasonable.

Best,

Jim


> targets:
> FileName	Rep	Treatment
> EF731D0.CEL	1	D
> EF731L0.CEL	1	L
> EF813D0.CEL	2	D
> EF813L0.CEL	2	L
> EF815D0.CEL	3	D
> EF815L0.CEL	3	L
>
> code:
> Rep = factor(targets$Rep) # this is like SibShip in the LIMMA example
> Treat = factor(targets$Treatment, levels=c("D","L"))
> design = model.matrix(~Rep+Treat)
> fit0 = lmFit(eset_t0, design)
> fit0 = eBayes(fit0)
> topTable(fit0, coef="TreatL", genelist=genetext, adjust="BH")
>
> Thanks in advance!
> Lauren
> _______________________________________________
> 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

--
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099



More information about the Bioconductor mailing list