# parametric scanning

parametric_scan <- function(Adata)
{

# Tomokazu Konishi, 27 Aug 2006�@�@v1.1, 19 Oct 2006
# as a supplement of an article entitled
# "Detection and Restoration of Hybridization Problems
# in Affymetrix GeneChip Data by Parametric Scanning"

### constants

ALL <- Adata@exprs
# subtraction of bkg
# for (i in 1:NumOfExp) ALL[,i] <- ALL[,i] -gamma[i]
samnam <- sampleNames(Adata)
NumOfExp <- length(samnam)
chip_row <- nrow(Adata); chip_col <- ncol(Adata)
row<-1:chip_row; col<-1:chip_col
NumOfData <- chip_row*chip_col
ideal_standard_hybridization <- array(NA, NumOfData)
dif <- array(NA, dim=c(NumOfData, NumOfExp))
logdif <- array(NA, dim=c(NumOfData,NumOfExp))
flag <- array(0, dim=c(chip_row,chip_col,NumOfExp))
sd_delta_z<-array(NA, dim=c(NumOfExp,1))
window_median <- array(0, dim=c(chip_row, chip_col, NumOfExp))

rl <- readline("Expected number of failures? ")
expect.failures <- as.numeric(rl)/2
save(expect.failures, file="expect_failures.rdata")
save(samnam, NumOfExp, chip_row, chip_col, row, col, NumOfData, file="constants.rdata")

## image generator function
save_image_png <- function(data, i, fname) {
w<-(chip_row+200); co<- (chip_col+300)
(fh <-chip_row/96); (fw <- chip_col/96)
png(filename = fname, width = w, height = co, pointsize = 20, bg="gray80", res = 96)
op <-par(pin = c(fh , fw))
image(row,col,data[row,col,i], zlim=c(-20,20), col=cm.colors(256), main=samnam[i])
par(op)
dev.off() }

## Student's normlization of measured data
for (i in 1:NumOfExp) {
medi <- median(ALL[,i], na.rm=T)
#cat(medi,"\t")
ALL[,i] <- log10(ALL[,i]/medi)
}

### ideal "standard hybridization"

for (j in 1:NumOfData) {
ideal_standard_hybridization[j] <- mean(ALL[j,], na.rm=T, trim=0.3) }
# the value
save(ideal_standard_hybridization, file="ideal_standard_hybridization.rdata")

### differences from the ideal standard
for (i in 1:NumOfExp) {
logdif[,i] <- ALL[,i] - ideal_standard_hybridization }

### robust z-normalization of differences

for (i in 1:NumOfExp) {
uq <- quantile(logdif[,i], 0.7, na.rm=T)
md <- quantile(logdif[,i], 0.5, na.rm=T)
dq <- quantile(logdif[,i] , 0.3, na.rm=T)
stdev <- (uq-dq)*(1/(qnorm(0.7)-qnorm(0.3)))
sd_delta_z[i, 1]<-stdev
logdif[,i] <- (logdif[,i]-md)/stdev
}
# sd_delta_z reflects the difference of each hybridization from the ideal standard
save(sd_delta_z, file="sd_delta_z.rdata")
# the value
save(logdif, file="logdif.rdata")

### pseudo image output
dim(logdif) <- c(chip_row, chip_col, NumOfExp)

for (i in 1:NumOfExp){
fname<- paste("image_", samnam[i],".png")
save_image_png(logdif, i, fname)
}

#### finding flags
### median of 5x5 window

cutoff<- qnorm(expect.failures/NumOfData, sd=0.31)*(-1)
# the value sd=0.31should be refined according to the mode of med_25_sd,
# which will be obtained in one of the checking sessions.
# the sd value will be reduced to 0.248 by progress of the GeneChip system.

for (i in 1:NumOfExp) {
for (c_row in 3:(chip_row-3)) {
for (c_col in 3:(chip_col-3)) {
b<-(logdif[(c_row-2):(c_row+2),(c_col-2):(c_col+2),i])
bna<-b[!is.na(b)]
if(length(bna)>(5*5*0.75)) {
bflag<-(median(bna, na.rm=T))
window_median[c_row, c_col, i]<-bflag
if (bflag > cutoff ) flag[(c_row-3):(c_row+3),(c_col-3):(c_col+3),i] <- NA
if (bflag < cutoff*(-1)) flag[(c_row-3):(c_row+3),(c_col-3):(c_col+3),i] <- NA
} } } }

dim(window_median)<-c(chip_row*chip_col, NumOfExp)
# to check distribution of the median of each windows' differences
save(window_median, file="window_median.rdata")

dim(flag) <- c(NumOfData,NumOfExp)
save(flag, file="flag.rdata")

#### reflection of the flags to the Adata
ALL <- Adata@exprs
ALL <- ALL+flag
Adata@exprs <- ALL
Adata<<-Adata
save(Adata, file="Adata.rdata")

### flag image output

dim(flag)<-c(chip_row,chip_col,NumOfExp)

for (i in 1:NumOfExp) {
fname<- paste("flag_", samnam[i],".png")
save_image_png (flag, i, fname)
}
#### end of the function
}
parametric_scan(Adata)