[BioC] use of duplicateCorrelation in Limma with agilent one-color arrays

Tue Feb 21 11:18:45 CET 2012

Thank you, Gordon, for your help. (coef obviously should be '1' not '-1')

Jabez

--- On Tue, 21/2/12, Gordon K Smyth <smyth at wehi.EDU.AU> wrote:

> From: Gordon K Smyth <smyth at wehi.EDU.AU>
> Subject: Re: use of duplicateCorrelation in Limma with agilent one-color arrays
> To: "Jabez Wilson" <jabezwuk at yahoo.co.uk>
> Cc: "Bioconductor mailing list" <bioconductor at r-project.org>
> Date: Tuesday, 21 February, 2012, 3:19
> Well, you could simplify to
> 
> Obj <- read.maimages(targets, source="agilent.median",
> green.only=TRUE)
> 
> and I recommend something more like offset=16 (any offset
> will give 
> positive results).  Also coef=-1 must surely give an
> empty result.
> 
> Best wishes
> Gordon
> 
> ---------------------------------------------
> Professor Gordon K Smyth,
> Bioinformatics Division,
> Walter and Eliza Hall Institute of Medical Research,
> 1G Royal Parade, Parkville, Vic 3052, Australia.
> smyth at wehi.edu.au
> http://www.wehi.edu.au
> http://www.statsci.org/smyth
> 
> On Mon, 20 Feb 2012, Jabez Wilson wrote:
> 
> > Dear Gordon, thank you for your comments. I am new to
> single color arrays, having used limma extensively with
> two-color in the past. I admit I have borrowed heavily with
> my coding from the mattick lab wiki(http://matticklab.com/index.php?title=Single_channel_analysis_of_Agilent_microarray_data_with_Limma),
> which you may have a view on.
> > I would be grateful if you could comment therefore on
> my amended script? (which with my data gives a
> corfit$consensus of ~0.4)
> > (also assuming the design matrix is correctly
> constructed)
> > Obj <- read.maimages(targets, source="agilent",
> green.only=TRUE, columns=list(G="gMedianSignal",
> Gb="gBGMedianSignal"))
> > Obj.corr <- backgroundCorrect (Obj,
> method="normexp", offset=1)
> > E <- normalizeBetweenArrays(Obj.corr,
> method="quantile")
> > #E.avg <= avereps(E, ID=E$genes$ProbeName)
> > then the analysis similar to before using E instead of
> E.avg
> > 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,13,13,14,14,15,15,16,16)
> > corfit <- duplicateCorrelation(E, design, ndups=1,
> block=biolrep)
> > fit <- lmFit(E, design, block=biolrep,
> cor=corfit$consensus)
> > cont.matrix <- makeContrasts("Treat1-WT",
> levels=design)
> > fit2 <- contrasts.fit(fit, cont.matrix)
> > fit2 <- eBayes(fit2)
> > topTable(fit2, coef=-1)
> >
> > Jabez
> >
> >
> >
> >
> > --- On Sun, 19/2/12, Gordon K Smyth <smyth at wehi.EDU.AU>
> wrote:
> >
> >> From: Gordon K Smyth <smyth at wehi.EDU.AU>
> >> Subject: use of duplicateCorrelation in Limma with
> agilent one-color arrays
> >> To: "Jabez Wilson" <jabezwuk at yahoo.co.uk>
> >> Cc: "Bioconductor mailing list" <bioconductor at r-project.org>
> >> Date: Sunday, 19 February, 2012, 10:20
> >> Dear Jabez,
> >>
> >> There a number of strange things about your code,
> which seem
> >> to be to do with trying to work around storing
> single color
> >> data in a 2-color data object.  Could you please
> use
> >> read.maimages() with green.only=TRUE, so that the
> data is
> >> read into a single color data object.  None of
> the
> >> work-around will then be necessary.
> >>
> >> You might also try computing the duplicate
> correlation
> >> without averaging replicate probes.
> >>
> >> Best wishes
> >> Gordon
> >>
> >>> Date: Fri, 17 Feb 2012 12:14:23 +0000 (GMT)
> >>> From: Jabez Wilson <jabezwuk at yahoo.co.uk>
> >>> To: bioconductor at r-project.org
> >>> Subject: [BioC] use of duplicateCorrelation in
> Limma
> >> with agilent
> >>>     one-color    arrays
> >>> Message-ID:
> >>>     <1329480863.2527.YahooMailClassic at web132503.mail.ird.yahoo.com>
> >>> Content-Type: text/plain; charset=iso-8859-1
> >>>
> >>> Hi, everyone, I am using limma to analyse an
> agilent
> >> one-color array experiment, and have run into
> difficulties
> >> with duplicateCorrelation.
> >>> My experiment is as follows: single color
> agilent
> >> arrays, 4 WT samples, and 3 samples of each of 4
> treatment
> >> (treatments 1-4). I also have technical replicates
> >> (replicated once) for each sample. There are
> therefore 32
> >> files. The targets file looks like this:
> >>> ?
> >>> ?? SampleNumber????????????????? FileName
> Condition
> >> Notes
> >>> 1???????????? 1? RH_02_1_77_Oct11_1_1.txt???
> Treat1??
> >> 1.1
> >>> 2???????????? 2? RH_02_1_77_Oct11_2_1.txt???
> Treat1??
> >> 1.2
> >>> 3???????????? 3? RH_04_1_77_Oct11_1_1.txt???
> Treat1??
> >> 2.1
> >>> 4???????????? 4? RH_07_1_77_Oct11_1_1.txt???
> Treat1??
> >> 2.2
> >>> 5???????????? 5? RH_04_1_77_Oct11_1_2.txt???
> Treat1??
> >> 3.1
> >>> 6???????????? 6? RH_07_1_77_Oct11_1_2.txt???
> Treat1??
> >> 3.2
> >>> 7???????????? 7? RH_04_1_77_Oct11_1_3.txt???
> Treat2??
> >> 4.1
> >>> 8???????????? 8? RH_07_1_77_Oct11_1_3.txt???
> Treat2??
> >> 4.2
> >>> 9???????????? 9? RH_04_1_77_Oct11_1_4.txt???
> Treat2??
> >> 5.1
> >>> 10?????????? 10? RH_07_1_77_Oct11_1_4.txt???
> Treat2??
> >> 5.2
> >>> 11?????????? 11? RH_04_1_77_Oct11_2_1.txt???
> Treat2??
> >> 6.1
> >>> 12?????????? 12? RH_07_1_77_Oct11_2_1.txt???
> Treat2??
> >> 6.2
> >>> 13?????????? 13 US0_05_1_77_Oct11_1_1.txt???
> Treat3??
> >> 7.1
> >>> 14?????????? 14? RH_01_1_77_Oct11_1_1.txt???
> Treat3??
> >> 7.2
> >>> 15?????????? 15 US0_05_1_77_Oct11_1_2.txt???
> Treat3??
> >> 8.1
> >>> 16?????????? 16? RH_01_1_77_Oct11_1_4.txt???
> Treat3??
> >> 8.2
> >>> 17?????????? 17? RH_02_1_77_Oct11_1_2.txt???
> Treat3??
> >> 9.1
> >>> 18?????????? 18? RH_02_1_77_Oct11_2_2.txt???
> Treat3??
> >> 9.2
> >>> 19?????????? 19 US0_05_1_77_Oct11_1_3.txt???
> Treat4?
> >> 10.1
> >>> 20?????????? 20? RH_01_1_77_Oct11_1_2.txt???
> Treat4?
> >> 10.2
> >>> 21?????????? 21 US0_05_1_77_Oct11_1_4.txt???
> Treat4?
> >> 11.1
> >>> 22?????????? 22? RH_01_1_77_Oct11_1_3.txt???
> Treat4?
> >> 11.2
> >>> 23?????????? 23? RH_04_1_77_Oct11_2_2.txt???
> Treat4?
> >> 12.1
> >>> 24?????????? 24? RH_07_1_77_Oct11_2_2.txt???
> Treat4?
> >> 12.2
> >>> 25?????????? 25
> US0_05_1_77_Oct11_2_1.txt??????? WT?
> >> 13.1
> >>> 26?????????? 26?
> RH_01_1_77_Oct11_2_1.txt??????? WT?
> >> 13.2
> >>> 27?????????? 27
> US0_05_1_77_Oct11_2_2.txt??????? WT?
> >> 14.1
> >>> 28?????????? 28?
> RH_01_1_77_Oct11_2_2.txt??????? WT?
> >> 14.2
> >>> 29?????????? 29
> US0_05_1_77_Oct11_2_3.txt??????? WT?
> >> 15.1
> >>> 30?????????? 30?
> RH_01_1_77_Oct11_2_3.txt??????? WT?
> >> 15.2
> >>> 31?????????? 31
> US0_05_1_77_Oct11_2_4.txt??????? WT?
> >> 16.1
> >>> 32?????????? 32?
> RH_01_1_77_Oct11_2_4.txt??????? WT?
> >> 16.2
> >>>
> >>> I run the following commands to process the
> data and
> >> create the design:
> >>> RG <- read.maimages(targets, columns =
> list(G =
> >> "gMedianSignal", Gb = "gBGMedianSignal", R =
> >> "gProcessedSignal",Rb = "gIsPosAndSignif"),
> annotation =
> >> c("Row",
> "Col","FeatureNum","ControlType","ProbeName"))
> >>> RG <- backgroundCorrect(RG,
> method="normexp",
> >> offset=1)
> >>> E <- normalizeBetweenArrays(RG,
> method="Aquantile")
> >>> E.avg <- avereps(E, ID=E$genes$ProbeName)
> >>> f <- factor(targets$Condition, levels =
> >> unique(targets$Condition))
> >>> design <- model.matrix(~0 + f)
> >>> colnames(design) <- levels(f)
> >>> ?
> >>> The problem arises when I do the
> duplicateCorrelation.
> >>> 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,13,13,14,14,15,15,16,16)
> >>> corfit <- duplicateCorrelation(E.avg,
> design,
> >> ndups=1,block=biolrep)
> >>> fit <- lmFit(E.avg$A, design,
> block=biolrep,
> >> cor=corfit$consensus.correlation)
> >>> contrast.matrix <-
> >>
> makeContrasts("Treat1-WT",levels=design)????????????????????????????????????????????????????????????fit2
> >> <- contrasts.fit(fit, contrast.matrix)
> >>> fit2 <- eBayes(fit2)
> >>> topTable(fit2, adjust="BH", coef="Treat1-WT",
> >> genelist=E.avg$genes, number=10)
> >>>
> >>> Whereas I would expect the
> corfit$consensus.correlation
> >> to be generally very positive, I get the value
> >> 0.01385223.?Does anyone have any suggestions? Any
> help would
> >> be appreciated
> >>> ?
> >>> Jabez
> >>>
> >>
> >>
> ______________________________________________________________________
> >> The information in this email is confidential and
> intended
> >> solely for the addressee.
> >> You must not disclose, forward, print or use it
> without the
> >> permission of the sender.
> >>
> ______________________________________________________________________
> >>
> >
> 
> ______________________________________________________________________
> The information in this email is confidential and inte...{{dropped:6}}