[BioC] Single Channel Analysis

[Gordon Smyth]

Dear Naomi,

At 06:20 AM 16/11/2005, Naomi Altman wrote:
>Hi Gordon,
>I just used the single channel analysis for a loop design.  I am wondering 
>if it is possible that there is a code error that swaps the roles of Cy3 
>and Cy5.  The results did not
>make sense when I fitted with as shown below.
>
>Cy5=c("S ","SB ","F ","St ","MB ","C ","P ","L ","C ",
>"MB ","St ","F ","SB ","S ","L ","P ")
>Cy3=c("MB ","C ","P ","L ","St ","S ","SB ","F ",
>"St","F ","P ","C ","MB ","L","SB","S")
>targets=data.frame(RG$targets,Cy3=Cy3,Cy5=Cy5)
>TargetArray=targetsA2C(targets)
>
>MAb <- normalizeBetweenArrays(MA,method="Aq")
>design <- model.matrix(~Target,data=TargetArray)
>corfit <- intraspotCorrelation(MAb,design)
>fit <-
>lmscFit(MAb,design,correlation=corfit$consensus)
>cont.matrix <- rbind(rep(0,7),diag(1,7))
>fit2 <- contrasts.fit(fit,cont.matrix)

Your contrasts compare all the other treatments vs "C", because "C" is 
alphabetically sorted first and becomes the first level of 
TargetArray$Target. Is that what you intended?
>But when I used
>targets=data.frame(RG$targets,Cy3=Cy5,Cy5=Cy3)
>
>the results make sense (we know what several of the genes do and they show 
>up in the right place this way) and the estimated error variance is cut in 
>half. I can tediously go through and compute the some of the tissue means 
>by hand to compare with the lmscFit output, but I thought you would know 
>where in the code to check this.
>
>Thanks,
>
>Naomi S. Altman                                814-865-3791 (voice)
>Associate Professor
>Dept. of Statistics                              814-863-7114 (fax)
>Penn State University                         814-865-1348 (Statistics)
>University Park, PA 16802-2111


The lmscFit() code has been well tested, and I like to think that such a 
drastic error as swapping Cy3 and Cy5 would have been eliminated right at 
the beginning. Here is a little analysis, taken from example(lmscFit) to 
show that lmscFit gives the same results as an ordinary log-ratio analysis, 
just with more significance.

Another similar example is at:
https://stat.ethz.ch/pipermail/bioconductor/2004-August/005607.html

Regards
Gordon

 > library(limma)
 > library(sma)
 > data(MouseArray)
 > RG <- backgroundCorrect(mouse.data, method = "half")
 > MA <- normalizeWithinArrays(RG, mouse.setup)
 > MA <- normalizeBetweenArrays(MA, method = "Aq")
 > targets <- data.frame(Cy3 = I(rep("Pool", 6)), Cy5 = I(c("WT", "WT", 
"WT", "KO", "KO", "KO")))
 >
 > # Log-ratio analysis
 > design <- cbind(WTvsPool=1,KOvsWT=c(0,0,0,1,1,1))
 > fit <- eBayes(lmFit(MA,design))
 >
 > # Separate channel analysis
 > targets.sc <- targetsA2C(targets)
 > targets.sc$Target <- factor(targets.sc$Target, levels = c("Pool", "WT", 
"KO"))
 > designsc <- model.matrix(~Target, data = targets.sc)
 > fitsc <- lmscFit(MA, designsc, correlation = 0.9)
 > cont.matrix <- cbind(KOvsWT = c(0, -1, 1))
 > fitsc2 <- contrasts.fit(fitsc, cont.matrix)
 > fitsc2 <- eBayes(fitsc2)
 >
 > topTable(fit,coef="KOvsWT")
               M        A          t     P.Value          B
2149 -2.9843122 12.29973 -17.730955 0.001779788 -0.1898774
5356 -1.8056831 12.80482  -9.588158 0.067261480 -0.6726402
2537 -1.0896738 13.59237  -8.988952 0.069578132 -0.7550846
1739 -1.1223897 13.79072  -6.746311 0.347698716 -1.2148770
1496 -0.9193937 12.24417  -6.330681 0.416827883 -1.3376946
1337 -1.1629891 11.63667  -5.995650 0.488202522 -1.4484874
540  -2.3767753 11.66888  -5.523154 0.692239877 -1.6252308
4139 -1.0183999 12.45938  -5.389478 0.694271543 -1.6800493
6319  0.8041058 10.95522   5.294233 0.694271543 -1.7205015
4941 -1.0259822 13.36684  -4.956517 0.837701802 -1.8738335
 > topTable(fitsc2)
               M         A          t      P.Value         B
2149 -2.9612419 12.299725 -17.654541 1.388509e-06 7.9167816
5356 -1.7778575 12.804818 -10.857724 2.453726e-04 5.7003388
2537 -1.0408212 13.592375  -7.813258 5.146981e-03 3.6076740
1739 -1.1411915 13.790725  -7.788565 5.146981e-03 3.5859380
1496 -0.8907386 12.244169  -6.858207 1.569256e-02 2.6981643
540  -2.3692707 11.668883  -6.204245 3.584067e-02 1.9879911
4139 -1.0215328 12.459384  -5.764137 6.258656e-02 1.4671351
4941 -1.0483388 13.366842  -5.022158 1.834206e-01 0.5085680
6319  0.9151930 10.955218   4.933243 1.834206e-01 0.3869647
5375  0.6655686  9.755678   4.924741 1.834206e-01 0.3752643