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

Mon Feb 20 16:03:34 CET 2012

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.
> ______________________________________________________________________
>