# Input files are tab-separated with this columns:
#
# Index	GridRow	GridCol	SpotRow	SpotCol	fgGreenMean	fgRedMean	bgGreenMean	bgGreenMedian	bgRedMean	bgRedMedian
#
#
# Usage example:
#
#
# m1 <- Microarray.init("c1SpotFile.dat");
#
# m1$scatterPlotRG();
# m1$scatterPlotMA();
# m1$drawRedArray();
# m1$drawGreenArray();
# m1$drawRGArray();
# m1$reportArray();
# m1$spotIntensity(channel = 2, log = TRUE, ratio = TRUE);
# m1$spotIntensity(channel = 1, log = TRUE, ratio = TRUE);

Microarray.init <- function(filename) {

	microrawdata <- read.table(filename, header = TRUE);
	colordepth=255;

	negfilter <- function(x) {
		if (x <= 0) NA  else x
	}

	spotIntensity <- function(channel) {
	      	color = if (channel == 1) "Green" else "Red";
        	intensity = microrawdata[[paste("fg",color,"Mean",sep="")]] - microrawdata[[paste("bg",color,"Mean",sep="")]];
        	intensity = mapply(negfilter, intensity);
        	intensity
	};

	getMatrixsub <- function(channel) {
		numgridrows = max(microrawdata[["GridRow"]]); numgridcols = max(microrawdata[["GridCol"]]); numspotrows = max(microrawdata[["SpotRow"]]); numspotcols = max(microrawdata[["SpotCol"]]);
		resmatrix = matrix(0, numspotrows * numgridrows,numspotcols * numgridcols);
        	intensities = spotIntensity(channel);
		for (i in 1: length(intensities)) {
			resmatrix[(microrawdata[["GridRow",i]] - 1) * numspotrows + microrawdata[["SpotRow",i]], (microrawdata[["GridCol",i]] - 1)*numspotcols + microrawdata[["SpotCol",i]]] = intensities[i];
		}
		resmatrix
	};

	greenmatrix <- getMatrixsub(1);
	redmatrix <- getMatrixsub(2);

	getMatrix <- function(channel = 1, logarithm = FALSE, ratio = FALSE) {
	    	res = if (ratio) greenmatrix/redmatrix else if (channel == 1) greenmatrix else redmatrix;
	    	if (logarithm) res = log2(res);
	    	res
	}

	drawArrayFromMatrix <- function(matrix, channel) {
	    	mycolors = if (channel == 1) rgb(0,0:colordepth,0,max=colordepth) else rgb(0:colordepth,0,0,max=colordepth);
	    	image(10 + (1:nrow(matrix)), 10 + (1:ncol(matrix)),matrix / max(matrix,na.rm=TRUE),col=mycolors, ann = FALSE, bg = "black");
	}

	drawRGArray <- function(greenmatrix, redmatrix) {
	    	mycolors = c();
		for (i in 0:colordepth) mycolors = c(mycolors,rgb(i,0:colordepth,0,max=colordepth));
		image(10 + (1:nrow(greenmatrix)), 10 + (1:ncol(greenmatrix)),floor((redmatrix/max(redmatrix,na.rm=TRUE))/(1/colordepth)) + greenmatrix/max(greenmatrix,na.rm=TRUE), zlim=c(0,(colordepth+1)), col = mycolors, ann = FALSE, bg = "black");
	};



	list (

		spotIntensity = function(channel = 1, logarithm = FALSE, ratio = FALSE) {
			res = if (ratio) spotIntensity(1)/spotIntensity(2) else spotIntensity(channel);
			if (logarithm) res = log2(res);
			res
		},

		getMatrix = getMatrix,
		drawGreenArray = function(logarithm = FALSE) {
                       drawArrayFromMatrix(getMatrix(1, logarithm = logarithm), 1);
                       title(main = paste("GreenArray : ", filename));
		},

		drawRedArray = function(logarithm = FALSE) {
                       drawArrayFromMatrix(getMatrix(1, logarithm = logarithm), 2);
                       title(main = paste("RedArray : ", filename));
		},

		drawRGArray = function(logarithm = FALSE) {
                       drawRGArray(getMatrix(1, logarithm = logarithm), getMatrix(2, logarithm = logarithm));
                       title(main = paste("RGArray : ", filename));
		},

		reportArray = function() {
			matr = log2(redmatrix) - log2(greenmatrix);
			cat("Median : ", median(matr, na.rm = TRUE), "\n");
			quant = quantile(matr, na.rm = TRUE);
			cat("Quartile : ", quant, "\n");
			cat("IQB :", quant[["75%"]] - quant[["25%"]], "\n");
			boxplot(matr);
			title(main = paste("Boxplot : ", filename))

		},

		scatterPlotRG = function() {
			plot(log2(redmatrix),log2(greenmatrix));
			title(main = paste("Scatterplot R/G : ", filename));
		},
		scatterPlotMA = function() {
			plot(((log2(redmatrix) + log2(greenmatrix))/2), (log2(redmatrix) - log2(greenmatrix)));
			title(main = paste("Scatterplot M/A : ", filename));
		}
	)
}