[BioC] limma: paired + multiple comparisons + technical replication?
Sarah Bonnin
Sarah.Bonnin at crg.eu
Thu Oct 3 18:07:18 CEST 2013
Dear Gordon,
Thank you for your answer.
The only thing is that I also want to compare all groups (1 to 1, all possibilities), and use the "makeContrasts" command - so as to automate better - and I'm not sure if my way of doing it is correct!
Here is my code, which actually works now, but I was wondering if it was really taking into account pairs as it should.
targets (a bit more informative than the previous one; TR is technical replicate)
FileName Groups Pairs
Sample1 Control 1
Sample1_TR Control 1
Sample2 Control 2
Sample3 Control 3
Sample4 Control 4
Sample5 Treat1 1
Sample5_TR Treat1 1
Sample6 Treat1 2
Sample7 Treat1 3
Sample8 Treat1 4
Sample9 Treat2 1
Sample9_TR Treat2 1
Sample10 Treat2 2
Sample11 Treat2 3
Sample12 Treat2 4
techrep <- c(1,1,2,3,4,5,5,6,7,8,9,9,10,11,12) # technical replication info
design <- model.matrix(~0+targets$Groups+targets$Pairs)
colnames(design) <- c(unique(targets$Groups),"Pairs")
corfit <- duplicateCorrelation(xx, design, ndups = 1, block = techrep)
fit <- lmFit(xx, design, block = techrep, cor = corfit$consensus)
contrast.matrix <- makeContrasts(contrasts=c("Treat1-Control", "Treat2-Control", "Treat2-Treat1"), levels=design)
fit <- contrasts.fit(fit, contrast.matrix)
fit <- eBayes(fit)
Thank you!
Sarah
-----Original Message-----
From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU]
Sent: Thursday, October 03, 2013 1:04 AM
To: Sarah Bonnin
Cc: Bioconductor mailing list
Subject: limma: paired + multiple comparisons + technical replication?
Dear Sarah,
Just put Pairs in the design matrix, for example
model.matrix(~Groups+Pairs)
and use duplicateCorrelation() for the biolreps.
Best wishes
Gordon
> Date: Tue, 1 Oct 2013 17:55:34 +0200
> From: Sarah Bonnin <Sarah.Bonnin at crg.eu>
> To: "bioconductor at r-project.org" <bioconductor at r-project.org>
> Subject: [BioC] limma: paired + multiple comparisons + technical
> replication ?
>
> Dear list,
>
> This question might be a bit redundant and I apologize for it, if it
> is, but I can't find a clear answer to what I'm trying to do.
>
> I am working on a set of 12 expression one-channel arrays.
>
> My target file is basically as follows:
> FileName Pairs Groups
> Sample1 1 Group1
> Sample2 1 Group1
> Sample3 1 Group2
> Sample4 1 Group2
> Sample5 1 Group3
> Sample6 1 Group3
> Sample7 2 Group1
> Sample8 2 Group2
> Sample9 2 Group3
> Sample10 3 Group1
> Sample11 3 Group2
> Sample12 3 Group3
>
> There are several parameters to take into account:
>
> - I want to produce all possible pairwise comparisons (Group3-Group2,
> Group2-Group1, Group3-Group1): "Groups" column
>
> - I want my design to take into account the paired samples: "Pairs"
> column
>
> - The last thing is that some samples are technical replicates
> (Sample1 with Sample2, Sample3 with Sample4, Sample5 with Sample6) and
> I would also like to take this into account.
>
> I've read the "8.7 Multi-level experiments" chapter of limma user guide, which guided me into combining paired data and multiple comparisons, in which case I would do:
>> design <- model.matrix(~0+factor(targets$Groups))
>> colnames(design) <- unique(targets$Groups) corfit <-
>> duplicateCorrelation(eset,design,block=targets$Pairs)
>> fit <-
>> lmFit(eset,design,block=targets$Pairs,correlation=corfit$consensus)
>
> In theory to take into account technical replicates I would use:
>> biolrep <- c(1,1,2,2,3,3,4,5,6,7,8,9) corfit <-
>> duplicateCorrelation(eset, block = biolrep) fit <- lmFit(eset, block
>> = biolrep, cor = corfit$consensus)
>
> But how can I combine all of this?
>
> Is there a way to somehow pass both paired and technical replication
> information into the duplicateCorrelation step? Or should I modify the
> design instead to take into account the paired design?
>
> It is getting quite confusing for me.
>
> Any help greatly appreciated!
>
> Thanks in advance!
>
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:6}}
More information about the Bioconductor
mailing list