[BioC] BioC] LIMMA: paired samples AND technical replicates
Gordon Smyth
smyth at wehi.EDU.AU
Fri Sep 21 02:39:41 CEST 2007
Charlotte,
It looks like it might work ok, but I don't have time to check it out properly.
Best wishes
Gordon
At 02:28 AM 21/09/2007, Charlotte Schjerling wrote:
>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