use of duplicateCorrelation in Limma with agilent one-color arrays
1
0
Entering edit mode
Jabez Wilson ▴ 150
@jabez-wilson-1839
Last seen 9.6 years ago
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
limma PROcess limma PROcess • 1.4k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 8 hours ago
WEHI, Melbourne, Australia

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

ADD COMMENT
0
Entering edit mode

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)

Jabe

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

Jabez

ADD REPLY

Login before adding your answer.

Traffic: 641 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6