# code for
# C. Kooperberg, T. G. Fazzio, J. J. Delrow, and T. Tskuiyama (2002)
# "Improved background correction for spotted DNA microarrays"
# Journal of Computational Biology, 9, 55-66.
# Copyright (C) 2001--2002 Charles Kooperberg
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# The text of the GNU General Public License, version 2, is available
# as http://www.gnu.org/copyleft or by writing to the Free Software
# Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
#
# 1) rewrite organize.data to fit your own data structure
# 2) run
# > orgdata <- organize.data(mydata)
# > bcorrecteddata <- bground(orgdata)
# 3) column 1 of bcorrecteddata has the background corrected data
#     on a log(2) scale (if log2=TRUE) or a regular scale (if log2=FALSE)
#     for the 635nm channel; column 2 has the same data for the 532nm channel
#     column 3 contains the quantity in equation (5) of Kooperberg et al. (2002)
#     for whichever channel this quantity was largest
# 4) note that the log2 option doesn't just take the log-base 2 of the data
#     from when log2=FALSE, as log2(E(X)) does not equal E(log2(X))
#
# Charles Kooperberg, 12/11/2002, clk@fhcrc.org
#


bground <- function(x,log2=TRUE)
{
        dx <- mtopfunct(x,mvarfactor(x),log2)[,c(4,11,18)]
        if(log2)dx[,1:2] <- dx[,1:2]/log(2)
        dx <- data.frame(dx)
        names(dx) <- c("x635","x532","warning")
        dx
}
organize.data <- function(datafile)
{
xx <- datafile[,c(1,2,3,9,11,12,14,18,20,21,23,34,35)]
#1#column 1 - block
#2#column 2 - column nr
#3#column 3 - row nr
#4#column 9 - foreground 635nm - median
#5#column 11 - foreground 635nm - SD
#6#column 12 - background 635nm - median
#7#column 14 - background 635nm - SD
#8#column 18 - foreground 532nm - median
#9#column 20 - foreground 532nm - SD
#10#column 21 - background 532nm - median
#11#column 23 - background 532nm - SD
#12#column 34 - forground pixels
#13#column 35 - background pixels
yy <- is.na(xx[,1])
for(i in 2:13) yy <- yy+is.na(xx[,i])
xx <- xx[yy==0,]
xx
}
mvarfactor <- function(x)
{
    xmm <- max(x[,1])
    b1 <- varaux1(x,6,xmm)
    b2 <- varaux1(x,10,xmm)
    c1 <- x[, 7]/sqrt(x[, 13])
    c2 <- x[, 11]/sqrt(x[, 13])
    m1 <- lm(b1 ~ c1 - 1, weights = 1/(c1 + 1))
    m2 <- lm(b2 ~ c2 - 1, weights = 1/(c2 + 1))
    c(m1$coef,m2$coef)
}
varaux1 <- function(x, j,xmm)
{
    uu <- varaux2(x, 1, j)
    for(i in 2:xmm) uu <- c(uu, varaux2(x, i, j))
    uu
}
varaux2 <- function(x, i, j)
{
    v1 <- x[x[, 1] == i, c(2, 3, j)]
    mm <- max(c(x[,2],x[,3]))+5
    v2 <- rep(-100, mm*mm)
    v3 <- 1 + v1[, 1] + v1[, 2] * mm
    v2[v3] <- v1[, 3]
    mm1 <- mm-1
    v4 <- matrix(v2, ncol = mm)
    v4a <- v4[c(-1, -mm), c(-1, -mm)]
    v4b <- v4[c(-1, -2), c(-1, -mm)]
    v4c <- v4[c(-1, -mm), c(-1, -2)]
    v4d <- v4[c(-mm1, -mm), c(-1, -mm)]
    v4e <- v4[c(-1, -mm), c(-mm1, -mm)]
    v4x <- cbind(as.vector(v4a), as.vector(v4b), as.vector(v4c), as.vector(v4d), as.vector(v4e))
    v4x[v4x > -100] <- v4x[v4x > -100] + 1
    v4x[v4x == -100] <- 0
    v5 <- apply(v4x, 1, sum)
    v6 <- apply(v4x != 0, 1, sum)
    v5 <- v5/v6
    v7 <- v4x - v5
    v7[v4x == 0] <- 0
    v8 <- apply(v7^2, 1, sum)/(v6 - 1)
    v9 <- matrix(0, ncol = mm, nrow = mm)
    v9[2:mm1, 2:mm1] <- v8
    v9[v4 <= 0] <- 0
    v10 <- as.vector(v9[v4 > 0])
    sqrt(v10)
}
mtopfunct <- function(x,mx,log2)
{
    y1 <- x[, c(4, 5, 12, 6, 7, 13)]
    y2 <- x[, c(8, 9, 12, 10, 11, 13)]
    y1[, c(2, 5)] <- y1[, c(2, 5)] * mx[1]
    y2[, c(2, 5)] <- y2[, c(2, 5)] * mx[2]
    t1 <- vtopaux1(y1,log2)
    t2 <- vtopaux1(y2,log2)
    x2 <- cbind(t1, t2)
    vv <- x2[, 7]
    vv[x2[, 14] > vv] <- x2[x2[, 14] > vv, 14]
    vv[vv < 0] <- 0
    x2 <- cbind(x2, x2[, 4] - x2[, 11], sqrt(x2[, 5]^2 + x2[, 12]^2),
        1 * (x2[ , 6] + x2[, 13] > 0), vv)
    x2
}
vtopaux1 <- function(x,log2)
{
    ww7 <- rep(0,32)
    yy7 <- rep(0,32)
    ww7[1 ]<- 0.00178328072169643; yy7[1 ]<- 0.99930504173577217;
    ww7[2 ]<- 0.00414703326056247; yy7[2 ]<- 0.99634011677195533;
    ww7[3 ]<- 0.00650445796897836; yy7[3 ]<- 0.99101337147674429;
    ww7[4 ]<- 0.00884675982636395; yy7[4 ]<- 0.98333625388462598;
    ww7[5 ]<- 0.01116813946013113; yy7[5 ]<- 0.97332682778991098;
    ww7[6 ]<- 0.01346304789671864; yy7[6 ]<- 0.96100879965205377;
    ww7[7 ]<- 0.01572603047602472; yy7[7 ]<- 0.94641137485840277;
    ww7[8 ]<- 0.01795171577569734; yy7[8 ]<- 0.92956917213193957;
    ww7[9 ]<- 0.02013482315353021; yy7[9 ]<- 0.91052213707850282;
    ww7[10]<- 0.02227017380838325; yy7[10]<- 0.88931544599511414;
    ww7[11]<- 0.02435270256871087; yy7[11]<- 0.86599939815409277;
    ww7[12]<- 0.02637746971505466; yy7[12]<- 0.84062929625258032;
    ww7[13]<- 0.02833967261425948; yy7[13]<- 0.81326531512279754;
    ww7[14]<- 0.03023465707240248; yy7[14]<- 0.78397235894334139;
    ww7[15]<- 0.03205792835485155; yy7[15]<- 0.75281990726053194;
    ww7[16]<- 0.03380516183714161; yy7[16]<- 0.71988185017161088;
    ww7[17]<- 0.03547221325688239; yy7[17]<- 0.68523631305423327;
    ww7[18]<- 0.03705512854024005; yy7[18]<- 0.64896547125465731;
    ww7[19]<- 0.03855015317861563; yy7[19]<- 0.61115535517239328;
    ww7[20]<- 0.03995374113272034; yy7[20]<- 0.57189564620263400;
    ww7[21]<- 0.04126256324262353; yy7[21]<- 0.53127946401989457;
    ww7[22]<- 0.04247351512365359; yy7[22]<- 0.48940314570705296;
    ww7[23]<- 0.04358372452932345; yy7[23]<- 0.44636601725346409;
    ww7[24]<- 0.04459055816375657; yy7[24]<- 0.40227015796399163;
    ww7[25]<- 0.04549162792741814; yy7[25]<- 0.35722015833766813;
    ww7[26]<- 0.04628479658131442; yy7[26]<- 0.31132287199021097;
    ww7[27]<- 0.04696818281621002; yy7[27]<- 0.26468716220876742;
    ww7[28]<- 0.04754016571483031; yy7[28]<- 0.21742364374000708;
    ww7[29]<- 0.04799938859645831; yy7[29]<- 0.16964442042399283;
    ww7[30]<- 0.04834476223480295; yy7[30]<- 0.12146281929612056;
    ww7[31]<- 0.04857546744150343; yy7[31]<- 0.07299312178779904;
    ww7[32]<- 0.04869095700913972; yy7[32]<- 0.02435029266342443;
    ww7 <- c(ww7,ww7)
    yy7 <- (c(yy7,-yy7))/2+0.5
    yy700 <- yy7
    for(i in 1:99)yy700 <- c(yy700,yy7+i)
    xx <- x[, 1]
    xy <- x[, 4]
    xs1 <- x[, 2]/sqrt(x[, 3])
    xs2 <- x[, 5]/sqrt(x[, 6])
    t1 <- rep(0, length(xx))
    t1 <- cbind(t1, t1, t1, t1, t1, t1, t1)
    for(i in 1:length(xx)) {
        t1[i, ] <- vtopaux2(xx[i], xy[i], xs1[i], xs2[i],ww7,yy700,log2)
        if(i/1000 == floor(i/1000))
            print(i)
    }
    t1
}
vtopaux2 <- function(x, y, sig1, sig2,ww7,yy7x,log2)
{
    upp <- 10 * sig1 + x
    nl <- 11
    wwh <- rep(ww7,nl)
    yyh <- (yy7x[1:(nl*length(ww7))])/nl
    b1 <- yyh * upp
    x1 <- vtopaux4(x, y, sig1, sig2)
    c1 <- (vtopaux3(b1, x, y, sig1, sig2) * upp)/x1
    c1w <- c1*wwh
    mwwh <- mean(wwh)
    m1 <- mean(c1w)/mwwh
    m1w <- m1*mwwh
    m2 <- mean(c1w * b1)/m1w
    tmp <- b1-m2
    m3 <- mean(c1w * tmp * tmp)/m1w
    if(log2){
        b2 <- log(b1)
    }
    else {
        b2 <- b1
    }
    m4 <- mean(c1w * b2)/(m1*mwwh)
    tmp <- b2-m4
    m5 <- sqrt(mean(c1w * tmp * tmp)/m1w)
    m6 <- 1 * ((abs(m1) - 1) > 0.05)
    m7 <- (y - x)/sqrt(sig1^2 + sig2^2)
    c(m1, m2, sqrt(m3), m4, m5, m6, m7)
}
vtopaux3 <- function(a, x, y, sig1, sig2)
{
    u <- x - a
    v <- y
    s <- sig1^2
    t <- sig2^2
    d10 <- 1/sqrt(s + t) * dnorm((u - v)/sqrt(s + t))
    d4 <- (u * t + v * s)/(s + t)
    d5 <- sqrt((s * t)/(s + t))
    d6 <- d10 * (1 - pnorm(0, d4, d5))
    d6
}
vtopaux4 <- function(x, y, sig1, sig2)
{
    upp <- 4 * sig2 + y
    b1 <- (((1:10000) - 0.5) * upp)/10000
    m1 <- mean(vtopaux5(b1, x, y, sig1, sig2)) * upp
    m1
}
vtopaux5 <- function(x, y1, y2, sig1, sig2)
{
    (1 - pnorm((x - y1)/sig1)) * dnorm((x - y2), sd = sig2)
}