[BioC] Another limma factorial that needs controlling and pairing

Karl Brand k.brand at erasmusmc.nl
Thu Jan 21 18:08:52 CET 2010


Dear BioC's,

I'd very much appreciate input on my design to achieve a sound analysis 
with limma.

Initially i thought i could get away with a two-factor design:

-treatments, 'S', 'E' & 'L'
-tissues, 'R' & 'C'

(targets file at bottom)

With 6 groups i adhered to the example in the limma user guide sec.8.7, 
thus-

targets <- readTargets("targets_no_time_no_animal.txt")
TP <- paste(targets$Tissue, targets$Pperiod, sep = ".")

TP <- factor(TP, levels = c("R.S", "C.S", "R.E", "C.E", "R.L", "C.L"))
design <- model.matrix(~0 + TP)
colnames(design) <- levels(TP)
fit <- lmFit(rma.pp, design)

cont.matrix <- makeContrasts(

    R.12t12vsR.6t18 = R.E - R.S,
    R.12t12vsR.18t6 = R.E - R.L,
    R.6t18vsR.18t6  = R.S - R.L,

    C.12t12vsC.6t18 = C.E - C.S,
    C.12t12vsC.18t6 = C.E - C.L,
    C.6t18vsC.18t6  = C.S - C.L,

    R.6t18vsC.6t18   = R.S - C.S,
    R.12t12vsC.12t12 = R.E - C.E,
    long.int.tis  = (R.E - R.S) - (C.E - C.S),
    short.int.tis = (R.E - R.L) - (C.E - C.L),

    levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit3 <- eBayes(fit2)
results <- (decideTests(fit3, method="global", adjust.method="BH", 
p.value=0.05))
etc.

This worked fine for my contrasts of interest until discussing with a 
statistician (thank you Renee) just what i was trying to get away with- 
within each of the 6 groups are 16 samples collected over a 24 hour 
period. These group time course's focus on another aspect of this 
experiment(*fyi see at bottom). What i'm trying to acheive now with 
limma is not concerned with 'time', merely identifying differentially 
expressed genes (thus not changing in time) between the 6 groups. 
However i understand this factor needs to be controlled for, not ignored 
as i did above.

And this is where my R & limma skills prove inadequate- controlling for 
the factor, Time.

Scouring the user guide and BioC list for examples of what i was 
attempting, i tried-

targets <- read.delim("RNA_Targets.txt")

    Tissue <- factor(targets$Tissue)
    Pperiod <- factor(targets$Pperiod)
    Time <- factor(targets$Time)

design <- model.matrix(~ Tissue + Pperiod + Time, data=targets)
colnames(design)

fit <- lmFit(rma.pp, design)
fit2 <- eBayes(fit)
topTable(fit2, coef=2)

_But_ i've failed to extract my contrasts of interest as detailed in the 
first design.

My primary appeal for help is then- how to obtain my contrast's of 
interest? Either from the second design, and/or how to setup an 
appropriate design which incorporates factor=time so that i can control 
for this?

Additionally, in attempting to really do the 'right thing' setting up 
this design, there is the 4th factor- pairing of the tissue's. As is 
evident form the targets "RNA_Targets.txt", each 'R' & 'C' pair come 
from the same animal. I expect incorporating this pairing into the model 
would add power. Any help/suggestions here would also be greatly 
appreciated.

Limma was built to easily handle all these parameters, but the limma/R 
expertise in this corner of the world are just too thin for me to take 
advantage of it w/o further help.

I look forward to any suggestions with thanks in advance, cheers,

Karl

"targets_no_time_no_animal.txt"- same as "RNA_Targets.txt" w/o columns 
'Time' & 'Animal'

"RNA_Targets.txt"-

FileName	Tissue	Pperiod	Time	Animal
01file.CEL	R	S	1	1
02file.CEL	C	S	1	1
03file.CEL	R	S	2	2
04file.CEL	C	S	2	2
05file.CEL	R	S	3	3
06file.CEL	C	S	3	3
07file.CEL	R	S	4	4
08file.CEL	C	S	4	4
09file.CEL	R	S	5	5
10file.CEL	C	S	5	5
11file.CEL	R	S	6	6
12file.CEL	C	S	6	6
13file.CEL	R	S	7	7
14file.CEL	C	S	7	7
15file.CEL	R	S	8	8
16file.CEL	C	S	8	8
17file.CEL	R	S	9	9
18file.CEL	C	S	9	9
19file.CEL	R	S	10	10
20file.CEL	C	S	10	10
21file.CEL	R	S	11	11
22file.CEL	C	S	11	11
23file.CEL	R	S	12	12
24file.CEL	C	S	12	12
25file.CEL	R	S	13	13
26file.CEL	C	S	13	13
27file.CEL	R	S	14	14
28file.CEL	C	S	14	14
29file.CEL	R	S	15	15
30file.CEL	C	S	15	15
31file.CEL	R	S	16	16
32file.CEL	C	S	16	16
33file.CEL	R	E	1	17
34file.CEL	C	E	1	17
35file.CEL	R	E	2	18
36file.CEL	C	E	2	18
37file.CEL	R	E	3	19
38file.CEL	C	E	3	19
39file.CEL	R	E	4	20
40file.CEL	C	E	4	20
41file.CEL	R	E	5	21
42file.CEL	C	E	5	21
43file.CEL	R	E	6	22
44file.CEL	C	E	6	22
45file.CEL	R	E	7	23
46file.CEL	C	E	7	23
47file.CEL	R	E	8	24
48file.CEL	C	E	8	24
49file.CEL	R	E	9	25
50file.CEL	C	E	9	25
51file.CEL	R	E	10	26
52file.CEL	C	E	10	26
53file.CEL	R	E	11	27
54file.CEL	C	E	11	27
55file.CEL	R	E	12	28
56file.CEL	C	E	12	28
57file.CEL	R	E	13	29
58file.CEL	C	E	13	29
59file.CEL	R	E	14	30
60file.CEL	C	E	14	30
61file.CEL	R	E	15	31
62file.CEL	C	E	15	31
63file.CEL	R	E	16	32
64file.CEL	C	E	16	32
65file.CEL	R	L	1	33
66file.CEL	C	L	1	33
67file.CEL	R	L	2	34
68file.CEL	C	L	2	34
69file.CEL	R	L	3	35
70file.CEL	C	L	3	35
71file.CEL	R	L	4	36
72file.CEL	C	L	4	36
73file.CEL	R	L	5	37
74file.CEL	C	L	5	37
75file.CEL	R	L	6	38
76file.CEL	C	L	6	38
77file.CEL	R	L	7	39
78file.CEL	C	L	7	39
79file.CEL	R	L	8	40
80file.CEL	C	L	8	40
81file.CEL	R	L	9	41
82file.CEL	C	L	9	41
83file.CEL	R	L	10	42
84file.CEL	C	L	10	42
85file.CEL	R	L	11	43
86file.CEL	C	L	11	43
87file.CEL	R	L	12	44
88file.CEL	C	L	12	44
89file.CEL	R	L	13	45
90file.CEL	C	L	13	45
91file.CEL	R	L	14	46
92file.CEL	C	L	14	46
93file.CEL	R	L	15	47
94file.CEL	C	L	15	47
95file.CEL	R	L	16	48
96file.CEL	C	L	16	48


*fyi- The group time course's, consisisting of 16 samples collected 
every 1.5 hours over 24 hours, attempt to identify 'circadian expression 
profiles' by F-testing the 16 samples for goodness-of-fit to a sinus 
curve. There is only one replicate per time point, relying instead on 
the density of samples, as opposed to replication. Clearly it would be 
ideal to incorporate this model into limma, as i understand is possible. 
But in failing to employ limma for 'routine' DE identification can't 
even dream of attempting such in the time available. Anyone who could do 
this would be a co-author on the paper. And *VERY* popular in the 
circadian field.

-- 
Karl Brand <k.brand-ampersand-erasmusmc.nl>
Department of Genetics
Erasmus MC
Dr Molewaterplein 50
3015 GE Rotterdam
lab +31 (0)10 704 3409 fax +31 (0)10 704 4743 mob +31 (0)642 777 268



More information about the Bioconductor mailing list