[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