Hi all,

                I might be misunderstanding something here, but I think I'm having trouble with Limma's duplicateCorrelation function.

I am analysing a set of custom spotted two color microarrays. Each comparison contains two arrays. Each array contains two copies of each spot. The design table for each comparison looks something like this:

>targets
                     SlideNumber             FileName   Cy3 Cy5       Date                 Name
yeast_wt_v_dd_rep1           823  s823_bwp17y+ddy.txt   wty ddy 17/12/2012   yeast_wt_v_dd_rep1

yeast_wt_v_dd_rep2           828     s828_wty+ddy.txt   wty ddy 18/01/2013   yeast_wt_v_dd_rep2


After background correction and normalisation, I run duplicateCorrelation before fitting the model:

>corfit <- duplicateCorrelation(MA_norm, design = c(1,1), ndups=2, spacing = 8324)

When I look at the results, I find that the consensus correlation is negative:


> corfit$consensus.correlation

[1] -0.4163954

Surely this is wrong? Examining the correlation for duplicate spots on each slide manually suggests a good correlation, particularly if one ignores flagged spots:


> good_spots <- MA_norm$weights[1:8324,] == 1 & MA_norm$weights[8325:16648,] ==1

> cor(MA_norm$M[1:8324,1][good_spots[,1]],MA_norm$M[8325:16648,1][good_spots[,1]])

[1] 0.8769604

> cor(MA_norm$M[1:8324,2][good_spots[,2]],MA_norm$M[8325:16648,2][good_spots[,2]])

[1] 0.858891

Even without removing the flagged spots, the correlation is still positive:

> cor(MA_norm$M[1:8324,1],MA_norm$M[8325:16648,1])

[1] 0.2669379

> cor(MA_norm$M[1:8324,2],MA_norm$M[8325:16648,2])

[1] 0.1843577

I know that this isn't the way that duplicateCorrelation calculates rho, but one would expect that a good correlation in this way might indicate at least a positive correlation from duplicateCorrelation? Or am I interpreting the consensus correlation incorrectly?

Here is the rest of my analysis script (I've left out some diagnostic plot generation):

>targets <- readTargets(row.names="Name")
>RG <- read.maimages(targets$FileName, source = "genepix", wt.fun=wtflags(weight=0, cutoff=-50))
>spottypes <- readSpotTypes()
>RG$genes$Status <- controlStatus(spottypes,RG)
>anno <- read.delim("anno.txt", comment.char="!", header = T)
>RG_background <- backgroundCorrect(RG,method = "normexp", offset =10)
>
>#flag out spots where intentity isn't much above background in both channels
>RG_background$weights[RG_background$R < 50 & RG_background$G < 50] = 0
>MA <- normalizeWithinArrays(RG_background)
>MA_norm <- MA[MA$genes$Status == "cDNA",]
>MA_norm$genes$Status <- NULL
>corfit <- duplicateCorrelation(MA_norm, design = c(1,1), ndups=2, spacing = 8324)
>fit <- lmFit(MA_norm, design = design, ndups = 2, correlation=corfit$consensus.correlation,spacing = 8324)
>fit<- eBayes(fit)

>sesssionInfo()

R version 2.15.1 (2012-06-22)

Platform: x86_64-pc-mingw32/x64 (64-bit)



locale:

[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252

[4] LC_NUMERIC=C                            LC_TIME=English_United Kingdom.1252



attached base packages:

[1] grid      stats     graphics  grDevices utils     datasets  methods   base



other attached packages:

 [1] dichromat_2.0-0   gplots_2.11.0     MASS_7.3-18       KernSmooth_2.23-7 caTools_1.14      gdata_2.12.0      gtools_2.7.0

 [8] reshape2_1.2.2    ggplot2_0.9.3     statmod_1.4.16    limma_3.14.4



loaded via a namespace (and not attached):

 [1] bitops_1.0-5       colorspace_1.2-1   digest_0.6.2       gtable_0.1.2       labeling_0.1       munsell_0.4        plyr_1.8

 [8] proto_0.3-10       RColorBrewer_1.0-5 scales_0.2.3       stringr_0.6.2      tools_2.15.1


Any advice appreciated

Ian


	[[alternative HTML version deleted]]

