Entering edit mode
Jabez Wilson
▴
150
@jabez-wilson-1839
Last seen 10.2 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
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)
then the analysis similar to before using E instead of E.avg
Jabe
Well, you could simplify to
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
Thank you, Gordon, for your help. (coef obviously should be '1' not '-1')
Jabez