library(affy)
library(affyPLM)
library(hgu133aprobe)
library(hgu133acdf)
library(hgu133a)
source("functions.R")
load("probes.rda")

attach(hgu133aprobe)
Index <- seq(along=Probe.Interrogation.Position)
tmp1 <- split(Probe.Interrogation.Position,Probe.Set.Name)
tmp2 <- split(Index,Probe.Set.Name)
Min <- lapply(tmp1,function(x) rep(min(x),length(x)))
Max <- lapply(tmp1,function(x) rep(max(x),length(x)))
start <- unlist(Min)[unlist(Index)]
end <- unlist(Max)[unlist(Index)]
pl <- (Probe.Interrogation.Position-start)/(end-start) ##proportion of location

rm(tmp1,tmp2,Min,Max)
gc()

tmp <- get("xy2i", "package:hgu133acdf")
pmIndex <- unlist(indexProbes(Data, "pm"))
subIndex <- match(tmp(x, y), pmIndex)


nData1 <- bg.correct.rma(Data[,which(pData(Data)[,4]=="one-round")])
nData1 <- normalize(nData1)

nData2 <- bg.correct.rma(Data[,which(pData(Data)[,4]=="two-round")])
nData2 <- normalize(nData2)

pms1 <- log2(pm(nData1)[subIndex,])
pms2 <- log2(pm(nData2)[subIndex,])

pd1 <- pData(nData1)
pd2 <- pData(nData2)

N11 <- sum(pd1[,3]=="testis/breast")
N12 <- sum(pd1[,3]=="testis/cervix")
N21 <- sum(pd2[,3]=="testis/breast")
N22 <- sum(pd2[,3]=="testis/cervix")
v1 <- ((N11-1)*rowSD(pms1[,pd1[,3]=="testis/breast"])^2+
       (N12-1)*rowSD(pms1[,pd1[,3]=="testis/cervix"])^2)/(N11+N12-2)
v2 <- ((N21-1)*rowSD(pms2[,pd2[,3]=="testis/breast"])^2+
       (N22-1)*rowSD(pms2[,pd2[,3]=="testis/cervix"])^2)/(N11+N12-2)

library(MASS)
Subset <- seq(1,length(v1),len=10000)
Y <- log((v1[order(pl)])[Subset])
X <- (sort(pl))[Subset]
#fit1 <- loess(Y~X,degree=1,span=1/3)
fit1 <- rlm(Y~X)
Y <- rowMeans(pms1)[order(pl)][Subset]
fit2 <- rlm(Y~X)
w1 <- rep(NA,length=length(probeNames(Data)))
w1[subIndex] <-  predict(fit2,newdata=data.frame(X=pl))/sqrt(exp(predict(fit1,newdata=data.frame(X=pl))))
w1[is.na(w1)] <- mean(w1,na.rm=TRUE)
w1 <- w1/mean(w1)
range(w1)

Subset <- seq(1,length(v2),len=10000)
Y <- log((v2[order(pl)])[Subset])
X <- (sort(pl))[Subset]
#fit1 <- loess(Y~X,degree=1,span=1/2)
fit1 <- rlm(Y~X)
Y <- rowMeans(pms2)[order(pl)][Subset]
fit2 <- rlm(Y~X)
w2 <- rep(NA,length=length(w1))
w2[subIndex] <- predict(fit2,newdata=data.frame(X=pl))/sqrt(exp(predict(fit1,newdata=data.frame(X=pl))))
w2[is.na(w2)] <- mean(w2,na.rm=TRUE)
w2 <- w2/mean(w2)
qqnorm(fit1$resid)
plot(X,fit1$fitted)
range(w2)



pns <- probeNames(nData2)
gns <- unique(pns)
eset <- matrix(NA,length(gns),nrow(pData(nData2)))
pms <- log2(pm(nData2))
for(i in seq(along=gns)){
  cat(i," ")
  Index <- which(pns==gns[i])
  PMS <- pms[Index,]
  probes <- as.factor(rep(1:nrow(PMS),ncol(PMS)))
  chips <- as.factor(rep(1:ncol(PMS),rep(nrow(PMS),ncol(PMS))))
  ws <- rep(w2[Index],ncol(PMS))
  fit1 <- rlm(as.vector(PMS)~-1+chips+probes,weight=ws)
  eset[i,] <- coef(fit1)[1:ncol(PMS)]
}
rownames(eset) <- gns
colnames(eset) <- rownames(pData(nData2))
WRMA2 <- new("exprSet",exprs=eset,phenoData=phenoData(nData2))
save(WRMA2,file="newexprs.rda")
             
pns <- probeNames(nData1)
gns <- unique(pns)
eset <- matrix(NA,length(gns),nrow(pData(nData1)))
pms <- log2(pm(nData1))
for(i in seq(along=gns)){
  cat(i," ")
  Index <- which(pns==gns[i])
  PMS <- pms[Index,]
  probes <- as.factor(rep(1:nrow(PMS),ncol(PMS)))
  chips <- as.factor(rep(1:ncol(PMS),rep(nrow(PMS),ncol(PMS))))
  ws <- rep(w1[Index],ncol(PMS))
  fit1 <- rlm(as.vector(PMS)~-1+chips+probes,weight=ws)
  eset[i,] <- coef(fit1)[1:ncol(PMS)]
}
rownames(eset) <- gns
colnames(eset) <- rownames(pData(nData1))
WRMA1 <- new("exprSet",exprs=eset,phenoData=phenoData(nData1))


save(WRMA1,WRMA2,w1,w2,file="newexprs.rda")