Hi,
I got an experimental design made of 3 replicates. Each replicate consists in a time course comparing the response of a cell line to a treatment "A" vs a placebo "B". Replicates are not carried out on the same day, but within each replicate, treatment and placebo time course are carried out in parallel, usually starting from the cell line culture.
Experimental questions are:
- How do the two treatments differ in their response over time?
- Do the two treatments differ overall?
I read the section "9.7 Multi-level Experiments" of limma user guide and Aaron's answer in LIMMA time series spline model.matrix design setup. I think there are two blocks (one for subjects, one for replicates), but I am not sure how to specify them. In the code below, I give my view of the analysis with only blocking for subjects.
How to setup the design and the analysis?
Any help is appreciated.
> cbind.data.frame(repl, trtm, timp, subj, g)[order(repl, trtm),]
repl trtm timp subj g
1 R1 A 1 AR1 A1
7 R1 A 2 AR1 A2
13 R1 A 3 AR1 A3
19 R1 A 4 AR1 A4
25 R1 A 5 AR1 A5
31 R1 A 6 AR1 A6
2 R1 B 1 BR1 B1
8 R1 B 2 BR1 B2
14 R1 B 3 BR1 B3
20 R1 B 4 BR1 B4
26 R1 B 5 BR1 B5
32 R1 B 6 BR1 B6
3 R2 A 1 AR2 A1
9 R2 A 2 AR2 A2
15 R2 A 3 AR2 A3
21 R2 A 4 AR2 A4
27 R2 A 5 AR2 A5
33 R2 A 6 AR2 A6
4 R2 B 1 BR2 B1
10 R2 B 2 BR2 B2
16 R2 B 3 BR2 B3
22 R2 B 4 BR2 B4
28 R2 B 5 BR2 B5
34 R2 B 6 BR2 B6
5 R3 A 1 AR3 A1
11 R3 A 2 AR3 A2
17 R3 A 3 AR3 A3
23 R3 A 4 AR3 A4
29 R3 A 5 AR3 A5
35 R3 A 6 AR3 A6
6 R3 B 1 BR3 B1
12 R3 B 2 BR3 B2
18 R3 B 3 BR3 B3
24 R3 B 4 BR3 B4
30 R3 B 5 BR3 B5
36 R3 B 6 BR3 B6
library(limma)
# design
trtm <- rep(c("A", "B"), 18)
timp <- rep(1:6, each=6)
repl <- rep(rep(c("R1", "R2", "R3"), each = 2), 6)
subj <- paste0(trtm, repl)
# g features points of interest to build contrasts
g <- factor(paste0(trtm, timp))
cbind.data.frame(trtm, timp, repl, subj, g)
# simulation + 1 experimental profile
med_sd = 0.04 # median sd observed
prfl = c( # observed means
A1=0.00, B1=0.20, A2=0.04, B2=0.60, A3=0.25, B3=0.95,
A4=0.10, B4=0.80, A5=0.15, B5=0.75, A6=0.15, B6=0.70)
set.seed(0)
meas = matrix(rnorm(250*length(g), sd = med_sd), ncol = length(g))
for (i in names(prfl)) {
meas[1, g==i] = meas[1, g==i] + prfl[i]
}
set.seed(1)
meas[-1,] = meas[-1,] + runif(nrow(meas)-1, 0, .3)
rownames(meas) = sprintf("G%03d", 1:nrow(meas))
dim(meas)
# graphical check
library(ggplot2)
ggplot(data.frame(trtm, timp, repl, subj, g, value=meas[1,]),
aes(timp, value, group=subj, col=repl)) + geom_point() + geom_line()
# limma design
design <- model.matrix(~ 0 + g)
colnames(design) <- levels(g)
# with blocking on subject in duplicateCorrelation().
corfit <- duplicateCorrelation(meas, design, block=subj)
corfit$consensus
# Fit the model
fit <- lmFit(meas, design, block=subj, correlation=corfit$consensus)
# How do the two treatments differ in their response over time?
con.1 <- makeContrasts(B1 - A1, levels=design)
fitcon.1 <- contrasts.fit(fit, con.1)
fitcon.1 <- eBayes(fitcon.1)
topTable(fitcon.1, number=5, adjust.method="BH")
# and so on...
# con.2 <- makeContrasts(B2 - A2, levels=design)
# con.3 <- makeContrasts(B3 - A3, levels=design)
# con.4 <- makeContrasts(B4 - A4, levels=design)
# con.5 <- makeContrasts(B5 - A5, levels=design)
# con.6 <- makeContrasts(B6 - A6, levels=design)
# Do the two treatments differ overall?
anova.con <- makeContrasts(
(B1 - A1), (B2 - A2), (B3 - A3),
(B4 - A4), (B5 - A5), (B6 - A6),
levels=design)
fitcon.anova <- contrasts.fit(fit, anova.con)
fitcon.anova <- eBayes(fitcon.anova)
topTable(fitcon.anova, number=5, adjust.method="BH")
# look at F test and p-value
# here we identify which tests to keep
rescon.anova = topTable(fitcon.anova, number=Inf, adjust.method="BH")
decideTests(fitcon.anova, method="nestedF")
Sincere thanks.