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

[Gordon K Smyth]

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 intend...{{dropped:4}}