[BioC] BioC] LIMMA: paired samples AND technical replicates
Charlotte Schjerling
araneus at mRNA.dk
Thu Sep 20 18:28:32 CEST 2007
Dear Gordon
I find my set-up a little more complicated than described in the
Limma manual so I have to ask whether the following is correct.
A text file with target information is read by:
targets <- readTargets()
> targets
SampleName Patient Treatment
1 AM92a A M.early
2 AM92b A M.early
3 AnM92a A nM.early
4 AnM92b A nM.early
5 AM01a A M.late
6 AM01b A M.late
7 AnM01a A nM.late
8 AnM01b A nM.late
9 BM92a B M.early
10 BM92b B M.early
11 BnM92a B nM.early
12 BnM92b B nM.early
13 BM02a B M.late
14 BM02b B M.late
15 BnM02a B nM.late
16 BnM02b B nM.late
17 DM88a C M.early
18 DM88b C M.early
19 DnM88a C nM.early
20 DnM88b C nM.early
21 DM97a C M.late
22 DM97b C M.late
23 DnM97a C nM.late
24 DnM97b C nM.late
Patient <- factor(targets$Patient)
Treat <- factor(targets$Treatment)
I do not want an intercept in my design matrix so I can either choose:
design <- model.matrix(~0+Patient+Treat)
> design
PatientA PatientB PatientC TreatM.late TreatnM.early TreatnM.late
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 1 0
4 1 0 0 0 1 0
5 1 0 0 1 0 0
6 1 0 0 1 0 0
7 1 0 0 0 0 1
8 1 0 0 0 0 1
9 0 1 0 0 0 0
10 0 1 0 0 0 0
11 0 1 0 0 1 0
12 0 1 0 0 1 0
13 0 1 0 1 0 0
14 0 1 0 1 0 0
15 0 1 0 0 0 1
16 0 1 0 0 0 1
17 0 0 1 0 0 0
18 0 0 1 0 0 0
19 0 0 1 0 1 0
20 0 0 1 0 1 0
21 0 0 1 1 0 0
22 0 0 1 1 0 0
23 0 0 1 0 0 1
24 0 0 1 0 0 1
attr(,"assign")
[1] 1 1 1 2 2 2
attr(,"contrasts")
attr(,"contrasts")$Patient
[1] "contr.treatment"
attr(,"contrasts")$Treat
[1] "contr.treatment"
or
design <- model.matrix(~0+Treat+Patient)
> design
TreatM.early TreatM.late TreatnM.early TreatnM.late PatientB PatientC
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 0 0 1 0 0 0
4 0 0 1 0 0 0
5 0 1 0 0 0 0
6 0 1 0 0 0 0
7 0 0 0 1 0 0
8 0 0 0 1 0 0
9 1 0 0 0 1 0
10 1 0 0 0 1 0
11 0 0 1 0 1 0
12 0 0 1 0 1 0
13 0 1 0 0 1 0
14 0 1 0 0 1 0
15 0 0 0 1 1 0
16 0 0 0 1 1 0
17 1 0 0 0 0 1
18 1 0 0 0 0 1
19 0 0 1 0 0 1
20 0 0 1 0 0 1
21 0 1 0 0 0 1
22 0 1 0 0 0 1
23 0 0 0 1 0 1
24 0 0 0 1 0 1
attr(,"assign")
[1] 1 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$Treat
[1] "contr.treatment"
attr(,"contrasts")$Patient
[1] "contr.treatment"
However, I think that only the latter can give me the contrasts that I want:
> contrasts
Contrasts
Levels TreatM.early - TreatnM.early TreatM.late -
TreatnM.late TreatM.late - TreatM.early TreatnM.late - TreatnM.early
TreatM.early 1
0 -1 0
TreatM.late 0
1 1 0
TreatnM.early -1
0 0 -1
TreatnM.late 0
-1 0 1
PatientB 0
0 0 0
PatientC 0
0 0 0
>
Will it be correct to use: design <- model.matrix(~0+Treat+Patient)?
And then I continue with:
biolrep <- c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12)
corfit <- duplicateCorrelation(eset, design, ndups=1, block=biolrep)
fit <- lmFit(eset, design, ndups=1, block=biolrep, cor=corfit$consensus)
contrasts <- makeContrasts(TreatM.early-TreatnM.early,
TreatM.late-TreatnM.late, TreatM.early-TreatM.late,
TreatnM.early-TreatnM.late, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
Will that give me the paired set-up so each of the 4 conditions are
treated as originating from the 3 patients (A, B, D)?
Thanks,
Charlotte
At 01:14 AM 9/7/2007, you wrote:
>Dear Charlotte,
>
>Paired samples are normally handled as described in Section 8.3
>"Paired Samples" of the limma User's Guide. There shouldn't be any
>conflict between this and your technical replicates.
>
>Best wishes
>Gordon
>
>>Date: Wed, 05 Sep 2007 17:54:42 +0200
>>From: Charlotte Schjerling <araneus at mRNA.dk>
>>Subject: [BioC] LIMMA: paired samples AND technical replicates
>>To: bioconductor at stat.math.ethz.ch
>>Message-ID: <20070905155442.ADC4E5C29 at ns1-int.rh.dk>
>>Content-Type: text/plain; charset=us-ascii; format=flowed
>>
>>I have tried to search the BioC archives but have not found posts
>>dealing with both paired samples and technical replicates in a setup
>>like mine.
>>
>>I have 3 patients A, B, and C from which I have 4 samples: M.early,
>>nM.early, M.late, and nM.late. Each sample has a technical
>>replicate resulting in 24 arrays (Affymetrix).
>>
>>biolrep Patient M.early nM.early M.late nM.late
>>1 A 1 0 0 0
>>1 A 1 0 0 0
>>2 A 0 1 0 0
>>2 A 0 1 0 0
>>3 A 0 0 1 0
>>3 A 0 0 1 0
>>4 A 0 0 0 1
>>4 A 0 0 0 1
>>5 B 1 0 0 0
>>5 B 1 0 0 0
>>6 B 0 1 0 0
>>6 B 0 1 0 0
>>7 B 0 0 1 0
>>7 B 0 0 1 0
>>8 B 0 0 0 1
>>8 B 0 0 0 1
>>9 C 1 0 0 0
>>9 C 1 0 0 0
>>10 C 0 1 0 0
>>10 C 0 1 0 0
>>11 C 0 0 1 0
>>11 C 0 0 1 0
>>12 C 0 0 0 1
>>12 C 0 0 0 1
>>
>>I would like to find the following 4 contrasts:
>>
>>1) M.early - nM.early
>>2) M.late - nM.late
>>3) M.early - M.late
>>4) nM.early - nM.late
>>
>>biolrep <- c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12)
>>corfit <- duplicateCorrelation(eset, design, ndups=1, block=biolrep)
>>
>>If I ignore the pairing part I have constructed the following design matrix:
>>
>> > design
>> M.early nM.early M.late nM.late
>>1 1 0 0 0
>>2 1 0 0 0
>>3 0 1 0 0
>>4 0 1 0 0
>>5 0 0 1 0
>>6 0 0 1 0
>>7 0 0 0 1
>>8 0 0 0 1
>>9 1 0 0 0
>>10 1 0 0 0
>>11 0 1 0 0
>>12 0 1 0 0
>>13 0 0 1 0
>>14 0 0 1 0
>>15 0 0 0 1
>>16 0 0 0 1
>>17 1 0 0 0
>>18 1 0 0 0
>>19 0 1 0 0
>>20 0 1 0 0
>>21 0 0 1 0
>>22 0 0 1 0
>>23 0 0 0 1
>>24 0 0 0 1
>>attr(,"assign")
>>[1] 1 1 1 1
>>attr(,"contrasts")
>>attr(,"contrasts")$f
>>[1] "contr.treatment"
>>
>>
>>fit <- lmFit(eset, design, ndups=1, block=biolrep, cor=corfit$consensus)
>>contrasts <- makeContrasts(M.early-nM.early, M.late-nM.late,
>>M.early-M.late, nM.early-nM.late, levels=design)
>>fit2 <- contrasts.fit(fit, contrasts)
>>fit2 <- eBayes(fit2)
>>
>>However, I would like to include the pairing by patient into the
>>analyses and that I cannot figure out how to do.
>>
>>Hope someone can help,
>>
>>Charlotte
More information about the Bioconductor
mailing list