###############################################################################
# In summary this script calculates the robust z score for one or more plates
# of data.
#
# Args:
#   input.data.frame contains one or more columns of plate data. The first column 
#   must contain the well annotation ie A01, A02 etc. Each subsequent column represents
#   a plate of data where the column name is the plate name.
#
#   plate.conf.data.frame is a file of the following format:
#	  well, type (syntax must be either "poscontr" or "negcontr" or "empty"), name (name of control)
#
# Returns:
#    robust z score normalized values per plate.
#
# Author: kate gould
###############################################################################

use strict;
use warnings;
use IO::Handle;
use File::Temp qw/ tempfile tempdir /;
my $tdir = tempdir( CLEANUP => 0 );

# check to make sure having correct input and output files
my $usage = "usage: robustZScoreNormalization.pl [TABULAR.in] [TABULAR.in] [TABULAR.out] \n";
die $usage unless @ARGV == 3;

#get the input arguments
my $plateData = $ARGV[0];
my $plateConfig =  $ARGV[1];
my $zNorm =  $ARGV[2];

#open the input files
open (INPUT1, "<", $plateData) || die("Could not open file $plateData \n");
open (INPUT2, "<", $plateConfig) || die("Could not open file $plateConfig \n");
open (OUTPUT1, ">", $zNorm) || die("Could not open file $zNorm \n");

#variable to store the name of the R script file
my $r_script;

# R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values 	
# construct an R script file and save it in a temp directory
chdir $tdir;
$r_script = "robustZScoreNormalization.r";

open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n"; 
print Rcmd "
	#options(show.error.messages = FALSE);
	header <- read.table(\"$plateData\", head=F, sep=\"\\t\", comment=\"\", nrows=1);
	input.data.frame <- read.table(\"$plateData\", head=F, sep=\"\\t\", skip=1);
	plate.conf.data.frame <- read.table(\"$plateConfig\", head=T, sep=\"\\t\", comment=\"\");

	# assumed second column is type ie negcontr or poscontr and third column is name of control
	excluded.wells <- unique(plate.conf.data.frame[,1])
	
	robust.z.score.normalise <- function(x, y) return ((x-median(y))/mad(y))
	robust.z.score.normalisation.data.frame <- data.frame(well=input.data.frame[,1], stringsAsFactors=FALSE)
	
	for (i in 2:length(colnames(input.data.frame))){
		replicate.sample.values <- input.data.frame[((!(input.data.frame[,1] %in% excluded.wells))&(!(is.na(input.data.frame[,i])))),][,i]
		robust.z.score.normalisation.data.frame <- cbind(robust.z.score.normalisation.data.frame, round(sapply(input.data.frame[,i], FUN=robust.z.score.normalise, replicate.sample.values), 2))
	}
	names(robust.z.score.normalisation.data.frame) <- names(input.data.frame)
	write.table(header, file=\"$zNorm\", quote=F, sep=\"\\t\", row.names=F, col.names=F);
	write.table(robust.z.score.normalisation.data.frame, file=\"$zNorm\", quote=F, sep=\"\\t\", row.names=F, col.names=F, append=T);
	#eof\n";

close Rcmd;	

system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");

#close the input and output files
close(OUTPUT1);
close(INPUT1);
close(INPUT2);