comparison mtls_analyze/heatmap.R @ 4:b465306d00ba draft default tip

Uploaded

4:b465306d00ba

GeneratePeakMatrix <- function(experiments, scores) {
  # Generates a score matrix for all mtls.
  #
  # Args:
  #   experiments: A list of underscore deliminated experiments for each mtl.
  #                There should never be a completely empty item in this list.
  #                for eg. experiments[1] = expA_expD_expB.
  #   scores: A list of underscore deliminated scores for each mtl, the length
  #           of these scores should be identicle to the length of the
  #           experiments. eg. scores[1] = 55_33_245.
  #
  # Returns:
  #   The peak score matrix for all mtls.
  experiments <- sapply(experiments, function(x) strsplit(x, split = "_"))
  scores <- sapply(scores, function(x) strsplit(x, split = "_"))
  unique.experiments <- unique(unlist(experiments))

  peaks=matrix(0,nr=length(experiments),nc=length(unique.experiments))
  colnames(peaks) <- unique.experiments

  for(i in 1:length(experiments)){
    for(j in 1:length(experiments[[i]])){
      peaks[i,experiments[[i]][j]]=as.numeric(scores[[i]][j])
    }
  }
  return(peaks)
##################################################################
}
GetChipData <- function(file.name,
                        proximity = "distal",
                        include.targetless = TRUE,
                        column.order = NA) {
  # Reads in, filters, and organizes mtls data
  #
  # Args:
  #   file.name: The path to the mtls file (including the file name).
  #   proximity: Either "distal" or "proximal", defines the gene target distance
  #              from the mtl. Default is "distal".
  #   include.targetless: If TRUE, includes mtls with no targets (applied after
  #                       the proximity filter); if not, mtls with no target
  #                       will be exclided. Default is TRUE.
  #   column.order: An optional vector of column names in the order in which
  #                 they will be presented. If this is left to default (NA) the
  #                 presented order of the chip columns will be the order in
  #                 which they are seen.
  #
  # Returns:
  #   An organized list of mtls data with the following elements:
  #     $peaks - a matrix of peak p-values
  #     $targets - a list of underscore deliminated gene targets for each mtl
  if(param$debug) {
    print("In GetChipData")
  }
  # Set Constants for the mtls file type:
  MTLNUM <- "mtl.id"
  CHROMOSOME <- "chr"
  EXPERIMENTS <- "expt"
  EXPERIMENTS.SORTED <- "expt.alphanum.sorted"
  START <- "start"
  END <- "end"
  SUMMIT <- "summit"
  SCORES <- "score"
  SCORE.MEAN <- "score.mean"
  SPANS <- "span.tfs"
  SPAN.TOTAL <- "span.l"
  PEAK.IDS <- "peak.ids"
  TARGETS <- "trgt"
#  PROXIMAL.TARGETS <- "trgt.prox"
#  DISTAL.TARGETS <- "trgt.dist"
#  PROXIMAL.TARGETS.TSS.DISTANCE <- "dtss.prox"
#  DISTAL.TARGETS.TSS.DISTANCE <- "dtss.dist"


  #Get chip data in from files:
  #TARGETS <- switch(proximity,
  #                   distal = DISTAL.TARGETS,
  #                   proximal = PROXIMAL.TARGETS,
  #                   stop(" Bad proximity argument supplied."))  # Default
cat("param$rna.files = ", param$rna.files, "\n")
if(!is.na(param$rna.files) && param$rna.files != "none") {
  keep.columns <- c(EXPERIMENTS, SCORES, TARGETS)
} else {
  keep.columns <- c(EXPERIMENTS, SCORES)
}
  file <- read.delim(file.name, header=T, as.is=T)[, keep.columns]
if(!is.na(param$rna.files) && param$rna.files != "none") {
  if (!include.targetless) {
    ix.has.trgt <- which(as.character(file[, TARGETS])!="")
    file <- file[ix.has.trgt, ]
  }
}
  chip = list()

  chip$peaks <- GeneratePeakMatrix(file[, EXPERIMENTS], file[, SCORES])
  if(!is.na(column.order)) {
  #If you specify an order, or a subset of experiments to include, this is where
  #that gets done.
  order <- unlist(strsplit(column.order, split="::"))
  chip$peaks <- chip$peaks[,order]
  }
  if(!is.na(param$rna.files) && param$rna.files != "none") {
  chip$targets <- file[ , TARGETS]
}
  
  return (chip)
}
GetRNADataFromOneFile <- function(file.name, fpkm = "avg") {
  # Reads in rna expression data from cufflinks
  #
  # Args:
  #   file.name: The path to the cufflinks file (including the file name).
  #   fpkm: What fpkm is chosen, options are hi, low and avg
  # Returns:
  #   An organized vector of rnaseq data with the following elements:
  #     names(data) - the genes from the file. eg names(data[1]) = JUND
  #     data - the fpkm score for that gene data[1] = 31 
  if(param$debug) {
    cat("In GetRNADataFromOneFile for ",file.name)
  }

 # Set Constants for the cufflinks file type:
  GENE <- "gene_id"
  BUNDLE <- "bundle_id"
  CHROMOSOME <- "chr"
  LEFT <- "left"
  RIGHT <- "right"
  FPKM_AVG <- "FPKM"
  FPKM_LOW <- "FPKM_conf_lo"
  FPKM_HIGH <- "FPKM_conf_hi"
  STATUS <- "status"

  
  #Get chip data in from files:
   FPKM <- switch(fpkm,
                  hi = FPKM_HIGH,
                  avg = FPKM_AVG,
                  low = FPKM_LOW,
                  stop(" Bad fpkm quality argument supplied."))  # Default
  
  keep.columns <- c(GENE, FPKM) 
  
  file <- read.delim(file.name, header=T, as.is=T)[, keep.columns]
  
  rna = vector(mode = "numeric", length = nrow(file))
  
  rna <- file[, FPKM]
  names(rna) <- file[ , GENE] 
  return (rna)
}

GetRNAData <- function(file.names, file.lables = file.names, fpkm = "avg") {
  # Reads in rna expression data from cufflinks
  #
  # Args:
  #   file.names: The list of paths to the cufflinks files (including the file name).
  #   fpkm: What fpkm is chosen, options are hi, low and avg
  # Returns:
  #   A matrix of rnaseq data with the following description:
  #     each col corresponds to an rnaseq run
  #     eacg row corresponds to a gene 
  if(param$debug) {
    print("In GetRNAData")
  }
  files <- list()
  for(i in 1:length(file.names)) {
  files[[i]] <-GetRNADataFromOneFile(file.names[i])
  }
  genes <- unique(names(unlist(files)))
  scores <- matrix(0,nrow=length(genes),ncol=length(file.names))
  rownames(scores) <- genes
  colnames(scores) <- file.names

  for(j in 1:length(file.names)) {
    scores[names(files[[j]]),j] <- files[[j]]
  }
print("# of cols for scores is ")
print(ncol(scores))
print ("file lables are ")
print (file.lables)
  colnames(scores) <- file.lables
  return(scores)
  #scores[genes,file.names] <- files[[]]
}

NormalizeRNA <- function(scores) {

#add psudocount
scores <- scores+1

numerator <- scores[,1:(ncol(scores)/2)]
denominator <- scores[,(ncol(scores)/2 + 1):ncol(scores)]

#new.scores <- scores[,-which(colnames(scores) == norm.exp)]
#norm.score <- scores[,which(colnames(scores) == norm.exp)]
#return(log2(new.scores/norm.score))

return(log2(numerator/denominator))
  
}

PrepareRNAforHeatmap <- function(new.scores.foldchange){
new.scores.split <- matrix(0, nr=nrow(new.scores.foldchange), nc=2*ncol(new.scores.foldchange))
is.even <- function(x){ x %% 2 == 0 }
corresponding.col <- function(x){ceiling(x/2)}
new.col.names <- vector(length = ncol(new.scores.split))
rownames(new.scores.split) <- rownames(new.scores.foldchange)
  for(i in 1:ncol(new.scores.split)) {
	if(is.even(i)) {sign <- "down"} else {sign <- "up"}
  new.col.names[i] <- paste(colnames(new.scores.foldchange)[corresponding.col(i)], 
					sign, sep =".")
	new.scores.split[,i] <- new.scores.foldchange[,corresponding.col(i)]
	if(is.even(i)) {
		new.scores.split[which(new.scores.split[,i]>0),i] <- 0
		new.scores.split[,i] <- -1*new.scores.split[,i]
	} else {
		new.scores.split[which(new.scores.split[,i]<0),i] <- 0
	}
	colnames(new.scores.split) <- new.col.names
}
return(new.scores.split)}

#ConvertExpressionToPval = function(expression, threshold = 1)
#{
#    print("In ConvertExpressionToPval")
#    #Generate Normal Stats
#    # expression.mean = apply(expression, 2, mean)
#    # expression.sd = apply(expression, 2, sd)
#    
#    ix.above.threshold = which(expression > threshold)
#    expression.mean = apply(expression, 2, function(i) mean(i[which(i > threshold)]))
#    expression.sd = apply(expression, 2, function(i)     sd(i[which(i > threshold)]))
#    

#    # CDF from zero to point
#    # expression.pval = matrix(,nr=nrow(expression),nc=ncol(expression),dimnames=dimnames(expression))
#    # ix.low = which(expression < expression.mean)
#    # ix.upper = which(expression >= expression.mean)
#    # expression.pval[ix.low] = (-10 * pnorm(expression[ix.low], mean = expression.mean, sd = expression.sd, log=T) *log10(exp(1)))
#    
#    expression.pval = (-10 * log10(exp(1)) * pnorm(expression, 
#                                                    mean = expression.mean, 
#                                                    sd = expression.sd, 
#                                                    log=T, 
#                                                    lower.tail=F
#    ) 
#    )
#    #squash those under threshold to zero
#    expression.pval[-ix.above.threshold] = 0

#    # correct for upper tail (ie from point to inf)
#    #Still need to do
#    return (expression.pval)
#}
MapExpressiontoChip2 = function(chip.targets, expression)
{
mymatrix = matrix(0, nr=length(chip.targets), nc=ncol(expression))
rownames(mymatrix) <- chip.targets
colnames(mymatrix) <- colnames(expression)


j <- 1 #expression counter
i <- 1 #mymatrix counter
expression.genes <- rownames(expression)
previous.expression.gene <- expression.genes[j]

# Progress bar:
total <- length(chip.targets)
# create progress bar
pb <- txtProgressBar(min = 0, max = total, style = 3)

for ( i in 1:length(chip.targets)) { 
  current.target <- chip.targets[i]
  setTxtProgressBar(pb, i)
  while (current.target > previous.expression.gene && j<=length(expression.genes)) {
    j <- j + 1
    previous.expression.gene <- expression.genes[j]
  }
  if (current.target == previous.expression.gene){
    mymatrix[i,] <- expression[j,]
  }
}

close(pb)

           print("Leaving MapRNAtoChip")
           return(mymatrix)
}


MapExpressiontoChip = function(chip.targets, expression)
{
    print("In MapRNAtoChip")
    # targets = unlist(strsplit(chip$targets[i], "_"))
    #   exp = matrix(rna$all.pval[targets, ])   
    #   return (apply(exp, 2, mean))
    
    #rna col 1 = th0 overexpressed  
    #rna col 2 = th17 overexpressed
    #rna col 3 = true zscores






    ####################################################################
    chip.targets.notnull.unique <- unique(chip.targets[which(chip.targets!="")])
    gene.intersect <- intersect(rownames(expression), chip.targets.notnull.unique)
    chip.targets.notnull.unique.with.data <- chip.targets.notnull.unique[which(chip.targets.notnull.unique %in% gene.intersect)]
    expression.useful.data <- expression[which(rownames(expression) %in% gene.intersect),]
    chip.targets.notnull.unique.with.data.sorted <- chip.targets.notnull.unique.with.data[order(chip.targets.notnull.unique.with.data)]
    expression.useful.data.sorted <- expression.useful.data[order(rownames(expression.useful.data)),]
    ####################################################################
    
    mymatrix = matrix(0, nr=length(chip.targets), nc=ncol(expression))
       rownames(mymatrix) <- chip.targets
       colnames(mymatrix) <- colnames(expression)
       head(chip.targets.notnull.unique.with.data.sorted)
       head(rownames(expression.useful.data.sorted))
       if(!identical(chip.targets.notnull.unique.with.data.sorted, rownames(expression.useful.data.sorted))){
        stop("We have a serious problem, chip is not alligned to expression")
      }
       for(i in 1:length(chip.targets.notnull.unique.with.data.sorted))
          {
            mtls.ix <- which(chip.targets == chip.targets.notnull.unique.with.data.sorted[i])
            for(j in 1:length(mtls.ix)){
                mymatrix[mtls.ix[j],] <- expression.useful.data.sorted[i,]
              }
          }
           print("Leaving MapRNAtoChip")
           return(mymatrix)
    
   #return( sapply(chip$targets, function(x) apply(rna$all.pval[unlist(strsplit(chip$targets[3], "_")), ], 2, mean)) )
}





GenerateKMOrder <- function(data, km) {
  # generates the order of clusters from high to low mean values
  #
  # Args:
  #   data: A matrix of scores that have already been clustered.
  #   km: The K-means object generated from running kmeans on data. Another
  #       method could be used so long as it supplies a (km)$cluser list. Must
  #       have the same length as the number of rows in data
  #
  # Returns:
  #   cluster.order: The order in which the clusters should be displayed.

    km.cluster = km$cluster
    clusters = unique(km.cluster)
    clusters.avg = numeric()
    for(i in clusters) {
        clusters.avg = c(clusters.avg, mean(data[which(km.cluster == i), ]))
    }
    if(param$debug) {
      print ("straight clusters")
      print (clusters)
      print ("straigth average")
      print (clusters.avg)
      print ("ordered clusters")
      print (clusters[order(clusters.avg)])
      print("ordered average")
      print (clusters.avg[order(clusters.avg)])
    }
    return(clusters[rev(order(clusters.avg))])
}

OrderMTL <- function(data, km, cluster.order) {
  # Orders a matrix of data according to a clustering algorithm
  #
  # Args:
  #   data: A matrix of scores that have already been clustered.
  #   km: The K-means object generated from running kmeans on data. Another
  #       method could be used so long as it supplies a (km)$cluser list. Must
  #       have the same length as the number of rows in data
  #   cluster.order: The order in which the clusters should be displayed. 
  #                  for eg. km.order = c(2, 3, 1) would result in cluster 2 
  #                  being on top, then cluster 3 and lastly cluster 1.
  #
  # Returns:
  #   a list that contains 3 objects:
  #     list$data: the ordered version of the data.
  #     list$color.vector: a list of colors that should be assigned to each row.
  #     list$start.row: the starting row of each cluster in data.
  number.clusters <- length(cluster.order)
  cluster.colors <- sample(rainbow(number.clusters))
  
  # Set up return objects
  sorted.data <- matrix(,nr=nrow(data), nc=ncol(data))
  colnames(sorted.data) <- colnames(data)
  cluster.color.vector = vector(length=length(km$cluster))
  cluster.start.row = numeric(number.clusters)
  cluster.start.row[1]=1
  
  for (i in 1:number.clusters)
  {
      current.cluster = cluster.order[i]
      ix = which(km$cluster == current.cluster)
      current.cluster.range <- cluster.start.row[i]:(cluster.start.row[i]+length(ix)-1)
      sorted.data[current.cluster.range, ] = data[ix, ]
      cluster.color.vector[current.cluster.range] = cluster.colors[i]
      cluster.start.row[i+1] = (cluster.start.row[i]+length(ix))  
  }
  ret.list = list()
  ret.list$data = sorted.data
  ret.list$color.vector = cluster.color.vector
  ret.list$cluster.index = cluster.start.row
  return(ret.list)
}    


CreateHeatMap <- function(data,
                         km,
                         cluster.order, 
                         document.name,
                         document.type = "png",
                         number.colors = 30) {
  # Generates a heatmap image based for a matrix based on a clustering 
  # algorithm.
  #
  # Args:
  #   data: A matrix of scores that have already been clustered, the column
  #         names of this matrix will become the column titels of the heatmap.
  #   km: The K-means object generated from running kmeans on data. Another
  #       method could be used so long as it supplies a (km)$cluser list.
  #   cluster.order: The order in which the clusters should be displayed. 
  #                  for eg. km.order = c(2, 3, 1) would result in cluster 2 
  #                  being on top, then cluster 3 and lastly cluster 1.
  #   document.name: A name for the produced file. there is no need to
  #                  supply the .png/.pdf in your argument.
  #   document.type: The type of file you want to produce. current options are
  #                  png and pdf. Default is pdf.
  #
  # Returns:
  #   Nothing to the script that calls it, however it creates an image at the 
  #   path specified.
  if(param$debug) {
    print("In CreateHeatMap")    
  }
  
  data.ordered <- OrderMTL(data, km, cluster.order)
    
	
	#Load Lib
	library(gplots)
	
	#Set Color gradient
	color.ramp = colorRampPalette(c("black",
	                                "darkblue",
	                                "blue",
	                                "yellow",
	                                "orange",
	                                "red"))(number.colors) #7
	
	if(document.type == "png") {
          png(paste(document.name,".png", sep = ""),
              #width=15360,#1920,
              #height=204080,#2560,
              #res=500,
              antialias="gray")
        }
        else if(document.type == "pdf") {
        pdf(paste(document.name,".pdf", sep = ""))
        }
        else if(document.type == "tiff") {
        tiff(paste(document.name,".tiff", sep = ""),
        res=800,
        pointsize=2,
        width=1920,
        height=1920)
        }
        else {
        bitmap(paste(document.name,".bmp", sep = ""),
        height = 5,
        width = 5,
        res = 500)
        }
#op <- par(mar = rep(0, 4))
	
	heatmap.2(
		data.ordered$data, #data.ordered$data[,sort(colnames(data.ordered$data))],
		rowsep = data.ordered$cluster.index[-1],
		sepwidth = c(0.5, ncol(data)/100), 
		dendrogram = "none",
		Rowv = F,
		Colv = F,
		trace = "none",
		labRow = F, #sapply(seq(1:length(data.ordered$cluster.index)), toString),
		labCol = colnames(data.ordered$data), #sort(colnames(data.ordered$data)), 
		RowSideColors = data.ordered$color.vector, 
		keysize=0.6,
		key=F,
		col = color.ramp,
		cexCol = 0.8,
		cexRow = 0.8)
  
	dev.off()  

}

CreateIndividualHeatMap <- function(data,
                         km,
                         cluster.order,
                         color.ramp = colorRampPalette(c("black",
                                  "darkblue",
                                  "blue",
                                  "yellow",
                                  "orange",
                                  "red"))(30)) {
  # Generates a heatmap image based for a matrix based on a clustering 
  # algorithm.
  #
  # Args:
  #   data: A matrix of scores that have already been clustered, the column
  #         names of this matrix will become the column titels of the heatmap.
  #   km: The K-means object generated from running kmeans on data. Another
  #       method could be used so long as it supplies a (km)$cluser list.
  #   cluster.order: The order in which the clusters should be displayed. 
  #                  for eg. km.order = c(2, 3, 1) would result in cluster 2 
  #                  being on top, then cluster 3 and lastly cluster 1.
  #   document.name: A name for the produced file. there is no need to
  #                  supply the .png/.pdf in your argument.
  #   document.type: The type of file you want to produce. current options are
  #                  png and pdf. Default is pdf.
  #
  # Returns:
  #   Nothing to the script that calls it, however it creates an image at the 
  #   path specified.
  if(param$debug) {
    print("In CreateHeatMap")    
  }
  
  data.ordered <- OrderMTL(data, km, cluster.order)
    
  
  #Load Lib
  library(gplots)
  

  

  heatmap.3(
    data.ordered$data, #data.ordered$data[,sort(colnames(data.ordered$data))],
    rowsep = data.ordered$cluster.index[-1],
    sepwidth = c(0.5, nrow(data.ordered$data)/100), 
    dendrogram = "none",
    Rowv = F,
    Colv = F,
    trace = "none",
    labRow = F, #sapply(seq(1:length(data.ordered$cluster.index)), toString),
    labCol = colnames(data.ordered$data), #sort(colnames(data.ordered$data)), 
    #sourcRowSideColors = data.ordered$color.vector, 
    keysize=0.6,
    key=F,
    col = color.ramp,
    cexCol = 0.8,
    cexRow = 0.8)
   

}

ReadCommadLineParameters <- function(argument.names, command.line.arguments, optional = F) {
  # Reads the parameters from the command line arguments.
  #
  # Args:
  #   argument.names: A list of expected argument names.
  #   command.line.arguments: The list of recieved command line arguments
  #   optional: Are the areguments optional, or are they required, default is required
  #
  # Returns:
  #   The arguments for argument.names. As strings In that order.

    if(length(grep("--version",command.line.arguments))) {
        cat("version",script.version,"\n")
        q()
    }
    # Split command line arguments
    args <- sapply(strsplit(command.line.arguments, " "),function(i) i)
    vals <- character(length(argument.names))
    # split cmd.line to key and value pairs
    for(i in 1:length(argument.names)) {
        ix <- grep(argument.names[i], args)
        if(length(ix)>1) {
            stop("arg ",
                 argument.names[i],
                 " used more than once.  Bailing out...\n",
                 PrintParamError())
        }
        else if (length(ix)==0 && !optional) {
            stop("could not find ",
                 argument.names[i],
                 ". Bailing out...\n",
                 PrintParamError())
        }
        else if (length(ix)==0 && optional) {
        vals[i] <- NA
        }
        else {
        vals[i] <- args[ix+1]
        }
    }
    return(vals)
}

PrintParamError <- function(){
  # Prints the usage of the function, shows users what arguments they can use
  #
  # Args:
  #   param: A list that contains all the paramaters.
  #
  # Returns:
  #   A modified version of the param list, with the default values loaded.
    cat("
DESCRIPTIION:
    heatmap.R takes a ...

INPUT:
    1.--mtls_file: path to mtls file.\n
    2.--cluster_file: the destination path for the output cluster file.\n
    3.--chip_experiment_order: The order of desired chip experiments (optional).\n
    4.--heatmap_file: path for output heatmap image (no extension).\n
    5.--heatmap_type: choice of image format, currently support png, pdf, tiff and bmp (optional)\n
    6.--expression_file: list of expression files to be included in analysis (optional).\n
    7.--expression_name: lables for the expression files (optional).\n
    8.--normalization_file: a list of files to be used for normalization,
    they can be the same file, however the number of expression nominated
    normalization files must match the number of expression files (optional)\n
    9.--n_clusters: number of clusters\n
    10.--filter_percentage: percentage of mtls that will be analysed. Eg: if
    we make filter_percentage 30, we will take the union of the top mtls in
    mean, non-zero mean and variance (optional).\n

EXAMPLE RUN:
    Rscript heatmap.R
        --mtls_file path/to/mtls.xls
        --cluster_file path/to/output/cluster
        --chip_experiment_order tf1::tf2::tf5::tf3
        --heatmap_file path/to/output/heatmap
        --heatmap_type png
        --expression_file path/to/exp1::path/to/exp2
        --expression_name myexp1::myexp2
        --normalization_file path/to/exp3::path/to/exp3
        --n_clusters 13
        --filter_percentage 100
        --include_targetless yes
        --number_bins 30

    ")
}

#LoadDefaultParams <- function(param) {
#  # Loads default paramaters for the heatmap application
#  #
#  # Args:
#  #   param: A list that contains all the previous paramaters.
#  #
#  # Returns:
#  #   A modified version of the param list, with the default values loaded.

#script.version=0.1
#param$debug = F

## RNA data:
#param$rna.files = ""
#param$rna.normalization = "none"

## Filter:
#param$filter.percentage <- 100

## Clustering:
#param$clustering.number.of.clusters <- 13

## Heatmap:
#param$heatmap.document.name <- "heatmap"
#param$heatmap.document.type <- "png"
##Cluster Groups:
#param$cluster.groups.document.name <- "clusters"

#return(param)

#}



LoadParams <- function(cmd.args, args.nms, n.start.numeric, optional = F) {
  # Loads user defined paramaters for the heatmap application
  #
  # Args:
  #   cmd.args: The command line arguments given.
  #   arg.nms: A list of possible command arguments.
  #   n.start.numeric: the first argument that should be numeric (alway put
  #                        these last).
  #   optional: This specifies if the params in cmd.args are optional or
  #             required.
  # Returns:
  #   A list of values assigned to each argument.


  vals <- ReadCommadLineParameters(argument.names = args.nms,
                                 command.line.arguments = cmd.args,
                                 optional = optional)

   #check if numeric params are indeed numeric
  if(!optional) {
  for(i in n.start.numeric:length(vals)){
    if(is.na(as.numeric(vals[i]))){
      stop("arg ",args.nms[i]," is not numeric.  Bailing out...\n",print.error())
    }
  }
}
  return (vals)

}

#ValidateParams <- function(params) {
#    return(T)
#}


##########

LoadDebugParams <- function(param) {

    cmd.args <- c(
      "--mtls_file data/test/mtls.xls",
      "--cluster_file data/test/cluster",
      "--heatmap_file data/test/heatmap",
      "--heatmap_type bmp",
      "--n_clusters 13",
      "--filter_percentage 100",
      "--expression_file /home/kieran/code/scorchR-heatmap/data/expression/rna1.tabular::/home/kieran/code/scorchR-heatmap/data/expression/rna2.tabular",
      "--expression_name batf.ko.0::batf.ko.17", 
      "--normalization_file mean",
      "--chip_experiment_order ac::bc::cc::dc::ec",
      "--include_targetless yes",
      "--number_bins 20"
    )

args.nms <- c(            "--mtls_file",        #1
                "--cluster_file",     #2
                "--chip_experiment_order",    #3
                "--heatmap_file",     #4
                "--heatmap_type",         #5
                "--expression_file",    #6
                "--expression_name",    #7
                "--normalization_file",    #8
                "--include_targetless",    #9
                 "--n_clusters",     #10
                "--filter_percentage", #11
                "--number_bins")    #12


# vals has the same order as args.nms
vals <- LoadParams(cmd.args, args.nms, n.start.numeric = 10, optional = F)


# ChIP data:
param$annotated.macs.file <- vals[1]
param$chip.order <- vals[3]
# RNA data:
param$rna.files = vals[6]
param$rna.names = vals[7]
param$rna.normalization.file = vals[8]
param$include.targetless = vals[9]

# Filter:
param$filter.percentage <- as.numeric(vals[11])

# Clustering:
param$number.bins <- as.numeric(vals[12])
param$clustering.number.of.clusters <- as.numeric(vals[10])

# Heatmap:
param$heatmap.document.name <- vals[4]
param$heatmap.document.type <- vals[5]
#Cluster Groups:
param$cluster.groups.document.name <- vals[2]

return(param)

}



LoadOptionalParams <- function(param) {

cmd.args <- commandArgs(trailingOnly = T)

args.nms <- c(  "--chip_experiment_order",  #1
                "--expression_file",        #2
                "--expression_name",        #3
                "--normalization_file",     #4
                "--heatmap_type",           #5
                "--include_targetless",     #6
                "--filter_percentage",      #7
                "--number_bins")            #8


# vals has the same order as args.nms
vals <- LoadParams(cmd.args, args.nms, n.start.numeric = 7, optional = T)

# ChIP data:
param$chip.order <- if(!is.na(vals[1])){vals[1]}else{NA}

# RNA data:
param$rna.files <- if(!is.na(vals[2])){vals[2]}else{"none"}
param$rna.names <- if(!is.na(vals[3])){vals[3]}else{"none"}
param$rna.normalization.file <- if(!is.na(vals[4])){vals[4]}else{"no"}
param$include.targetless <- if(!is.na(vals[6])){vals[6]}else{"yes"}
# Filter:
param$filter.percentage <- if(!is.na(vals[7])){as.numeric(vals[7])}else{100}
param$number.bins <- if(!is.na(vals[8])){as.numeric(vals[8])}else{30}
# Heatmap file (output)
param$heatmap.document.type <- if(!is.na(vals[5])){vals[5]}else{"none"}

return(param)
}


LoadReqiredParams <- function(param){

cmd.args <- commandArgs(trailingOnly = T)

args.nms <- c(  "--mtls_file",          #1
                "--cluster_file",       #2
                "--heatmap_file",       #3
                 "--n_clusters")        #4


# vals has the same order as args.nms
vals <- LoadParams(cmd.args, args.nms, n.start.numeric = 4, optional = F)

# ChIP data:
param$annotated.macs.file <- vals[1]

# Clustering
param$clustering.number.of.clusters <- as.numeric(vals[4])

# Cluster file (output):
param$cluster.groups.document.name <- vals[2]


# Heatmap file (output):
param$heatmap.document.name <- vals[3]

return(param)

}

###########
# here we output the fasta file of the targets of each kmeans cluster
CreateClusterGroups <- function(trgts,k.ix,f.nm="output/clust.to.trgts", km.order){
  f.nm = paste(f.nm, ".fasta", sep="")
    clusters = km.order
    for(i in 1:length(clusters)){
        v=trgts[which(k.ix==clusters[i])]
        v.split = unlist(sapply(v,strsplit, "_"))
        if(i == 1){
            cat(sep="",file=f.nm,">cluster_",i,"\n")#clusters[i],"\n")    
        } else {
            cat(sep="",file=f.nm,">cluster_",i,"\n",append=TRUE)#clusters[i],"\n",append=TRUE)
        }
        cat(sep="\n",file=f.nm,v.split,append=TRUE)
    }
}


PrintClusters <- function(trgts,k.ix,f.nm="output/clust.to.trgts", km.order){
	f.nm = paste(f.nm, ".tsv", sep="")
	cat(sep="\t",file=f.nm,"row_number/target","cluster","\n")
	trgts[which(trgts=="")] <- "no_target"
    for(i in 1:length(trgts)){
		cat(sep="",file=f.nm,trgts[i],"\t","cluster_",k.ix[i],"\n",append=TRUE)
			}
#  	  if(i == 1){
#  	    cat(sep="",file=f.nm,"cluster_",i,"\n")#clusters[i],"\n")    
#  	  } else {
#  	    cat(sep="",file=f.nm,">cluster_",i,"\n",append=TRUE)#clusters[i],"\n",append=TRUE)
#  	  }
      #cat(sep="\n",file=f.nm,v.split,append=TRUE)
  	}
GetTopRowsFromMatrix = function(mtrx, percentage = 10)
{
  if (param$debug) {
    print("In GetTopRowsFromMatrix")
	}
	#Store the stats for the mtrx
	stats = list()
	stats$mean=apply(mtrx,1,mean)
	stats$sd=apply(mtrx,1,sd)
	stats$nonzero.mean=apply(mtrx, 1, function(x) mean(x[which(x != 0)]))
    
	
	#Store the indexes for the mtrx
	index = list()
	index$mean = sort(stats$mean, decreasing=T, index.return=T)$ix
	index$sd   = sort(stats$sd, decreasing=T, index.return=T)$ix
	index$nonzero.mean = sort(stats$nonzero.mean, decreasing=T, index.return=T)$ix
	
	#Calculate how many data points we want to take
	cutoff = floor(length(mtrx[ ,1])*(percentage/100))
	
	#Extract Indexes and return to caller
	index$union = union(index$mean[1:cutoff], union(index$sd[1:cutoff], index$nonzero.mean[1:cutoff]))
	
	return(index)
}

#GetUniq = function(targets.cluster)
#{
#    #problem with this function is that it agrigates the list handed to it. after this point you cant maintain order
#    
#    return(unique(unlist(lapply(targets.cluster, function(i) strsplit(i, "_")))))
#    #return(unlist(lapply(targets.cluster, function(i) strsplit(i, "_"))))
#    
#}

bindata = function(d,qunts=seq(.4,.9,.1))
{
    tmp=matrix(0,nr=nrow(d),nc=ncol(d),dimnames=dimnames(d))
    for(j in 1:ncol(d))
    {
        bins=quantile(d[,j],qunts)
            for(i in 1:length(bins))
            {
                tmp[which(d[,j]>bins[i]),j] = i
         }
    }
    return(tmp)
}

bindata.non.zero = function(d,qunts=seq(0,0.9,0.1))
{
    tmp=matrix(0,nr=nrow(d),nc=ncol(d))
    for(j in 1:ncol(d))
    {
        ix.non.zero=which(d[,j]!=0)
        bins=quantile(d[ix.non.zero,j],qunts)
            for(i in 1:length(bins))
            {
                tmp[which(d[,j]>bins[i]),j] = i
         }
    }
    
    return(tmp)
}
bindata.matrix = function(d,qunts=seq(0,0.9,0.1))
{
    tmp=matrix(0,nr=nrow(d),nc=ncol(d),dimnames=dimnames(d))
    #ix.non.zero=which(d!=0)
    bins=quantile(d[],qunts)
    for(i in 1:length(bins))
    {
        tmp[which(d>bins[i])] = i
    }
    return(tmp)
}
bindata.non.zero.matrix = function(d,qunts=seq(0,0.9,0.1))
{
    tmp=matrix(0,nr=nrow(d),nc=ncol(d),dimnames=dimnames(d))
    ix.non.zero=which(d!=0)
    bins=quantile(d[ix.non.zero],qunts)
    for(i in 1:length(bins))
    {
        tmp[which(d>bins[i])] = i
    }
    return(tmp)
}
heatmap.3 <- function (x, Rowv = TRUE, Colv = if (symm) "Rowv" else TRUE, 
    distfun = dist, hclustfun = hclust, dendrogram = c("both", 
        "row", "column", "none"), symm = FALSE, scale = c("none", 
        "row", "column"), na.rm = TRUE, revC = identical(Colv, 
        "Rowv"), add.expr, breaks, symbreaks = min(x < 0, na.rm = TRUE) || 
        scale != "none", col = "heat.colors", colsep, rowsep, 
    sepcolor = "white", sepwidth = c(0.05, 0.05), cellnote, notecex = 1, 
    notecol = "cyan", na.color = par("bg"), trace = c("column", 
        "row", "both", "none"), tracecol = "cyan", hline = median(breaks), 
    vline = median(breaks), linecol = tracecol, margins = c(5, 
        5), ColSideColors, RowSideColors, cexRow = 0.2 + 1/log10(nr), 
    cexCol = 0.2 + 1/log10(nc), labRow = NULL, labCol = NULL, 
    key = TRUE, keysize = 1.5, density.info = c("histogram", 
        "density", "none"), denscol = tracecol, symkey = min(x < 
        0, na.rm = TRUE) || symbreaks, densadj = 0.25, main = NULL, 
    xlab = NULL, ylab = NULL, lmat = NULL, lhei = NULL, lwid = NULL, 
    ...) 
{
    scale01 <- function(x, low = min(x), high = max(x)) {
        x <- (x - low)/(high - low)
        x
    }
    retval <- list()
    scale <- if (symm && missing(scale)) 
        "none"
    else match.arg(scale)
    dendrogram <- match.arg(dendrogram)
    trace <- match.arg(trace)
    density.info <- match.arg(density.info)
    if (length(col) == 1 && is.character(col)) 
        col <- get(col, mode = "function")
    if (!missing(breaks) && (scale != "none")) 
        warning("Using scale=\"row\" or scale=\"column\" when breaks are", 
            "specified can produce unpredictable results.", "Please consider using only one or the other.")
    if (is.null(Rowv) || is.na(Rowv)) 
        Rowv <- FALSE
    if (is.null(Colv) || is.na(Colv)) 
        Colv <- FALSE
    else if (Colv == "Rowv" && !isTRUE(Rowv)) 
        Colv <- FALSE
    if (length(di <- dim(x)) != 2 || !is.numeric(x)) 
        stop("`x' must be a numeric matrix")
    nr <- di[1]
    nc <- di[2]
    if (nr <= 1 || nc <= 1) 
        stop("`x' must have at least 2 rows and 2 columns")
    if (!is.numeric(margins) || length(margins) != 2) 
        stop("`margins' must be a numeric vector of length 2")
    if (missing(cellnote)) 
        cellnote <- matrix("", ncol = ncol(x), nrow = nrow(x))
    if (!inherits(Rowv, "dendrogram")) {
        if (((!isTRUE(Rowv)) || (is.null(Rowv))) && (dendrogram %in% 
            c("both", "row"))) {
            if (is.logical(Colv) && (Colv)) 
                dendrogram <- "column"
            else dedrogram <- "none"
            warning("Discrepancy: Rowv is FALSE, while dendrogram is `", 
                dendrogram, "'. Omitting row dendogram.")
        }
    }
    if (!inherits(Colv, "dendrogram")) {
        if (((!isTRUE(Colv)) || (is.null(Colv))) && (dendrogram %in% 
            c("both", "column"))) {
            if (is.logical(Rowv) && (Rowv)) 
                dendrogram <- "row"
            else dendrogram <- "none"
            warning("Discrepancy: Colv is FALSE, while dendrogram is `", 
                dendrogram, "'. Omitting column dendogram.")
        }
    }
    if (inherits(Rowv, "dendrogram")) {
        ddr <- Rowv
        rowInd <- order.dendrogram(ddr)
    }
    else if (is.integer(Rowv)) {
        hcr <- hclustfun(distfun(x))
        ddr <- as.dendrogram(hcr)
        ddr <- reorder(ddr, Rowv)
        rowInd <- order.dendrogram(ddr)
        if (nr != length(rowInd)) 
            stop("row dendrogram ordering gave index of wrong length")
    }
    else if (isTRUE(Rowv)) {
        Rowv <- rowMeans(x, na.rm = na.rm)
        hcr <- hclustfun(distfun(x))
        ddr <- as.dendrogram(hcr)
        ddr <- reorder(ddr, Rowv)
        rowInd <- order.dendrogram(ddr)
        if (nr != length(rowInd)) 
            stop("row dendrogram ordering gave index of wrong length")
    }
    else {
        rowInd <- nr:1
    }
    if (inherits(Colv, "dendrogram")) {
        ddc <- Colv
        colInd <- order.dendrogram(ddc)
    }
    else if (identical(Colv, "Rowv")) {
        if (nr != nc) 
            stop("Colv = \"Rowv\" but nrow(x) != ncol(x)")
        if (exists("ddr")) {
            ddc <- ddr
            colInd <- order.dendrogram(ddc)
        }
        else colInd <- rowInd
    }
    else if (is.integer(Colv)) {
        hcc <- hclustfun(distfun(if (symm) 
            x
        else t(x)))
        ddc <- as.dendrogram(hcc)
        ddc <- reorder(ddc, Colv)
        colInd <- order.dendrogram(ddc)
        if (nc != length(colInd)) 
            stop("column dendrogram ordering gave index of wrong length")
    }
    else if (isTRUE(Colv)) {
        Colv <- colMeans(x, na.rm = na.rm)
        hcc <- hclustfun(distfun(if (symm) 
            x
        else t(x)))
        ddc <- as.dendrogram(hcc)
        ddc <- reorder(ddc, Colv)
        colInd <- order.dendrogram(ddc)
        if (nc != length(colInd)) 
            stop("column dendrogram ordering gave index of wrong length")
    }
    else {
        colInd <- 1:nc
    }
    retval$rowInd <- rowInd
    retval$colInd <- colInd
    retval$call <- match.call()
    x <- x[rowInd, colInd]
    x.unscaled <- x
    cellnote <- cellnote[rowInd, colInd]
    if (is.null(labRow)) 
        labRow <- if (is.null(rownames(x))) 
            (1:nr)[rowInd]
        else rownames(x)
    else labRow <- labRow[rowInd]
    if (is.null(labCol)) 
        labCol <- if (is.null(colnames(x))) 
            (1:nc)[colInd]
        else colnames(x)
    else labCol <- labCol[colInd]
    if (scale == "row") {
        retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
        x <- sweep(x, 1, rm)
        retval$rowSDs <- sx <- apply(x, 1, sd, na.rm = na.rm)
        x <- sweep(x, 1, sx, "/")
    }
    else if (scale == "column") {
        retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
        x <- sweep(x, 2, rm)
        retval$colSDs <- sx <- apply(x, 2, sd, na.rm = na.rm)
        x <- sweep(x, 2, sx, "/")
    }
    if (missing(breaks) || is.null(breaks) || length(breaks) < 
        1) {
        if (missing(col) || is.function(col)) 
            breaks <- 16
        else breaks <- length(col) + 1
    }
    if (length(breaks) == 1) {
        if (!symbreaks) 
            breaks <- seq(min(x, na.rm = na.rm), max(x, na.rm = na.rm), 
                length = breaks)
        else {
            extreme <- max(abs(x), na.rm = TRUE)
            breaks <- seq(-extreme, extreme, length = breaks)
        }
    }
    nbr <- length(breaks)
    ncol <- length(breaks) - 1
    if (class(col) == "function") 
        col <- col(ncol)
    min.breaks <- min(breaks)
    max.breaks <- max(breaks)
    x[x < min.breaks] <- min.breaks
    x[x > max.breaks] <- max.breaks
    # if (missing(lhei) || is.null(lhei)) 
    #     lhei <- c(keysize, 4)
    # if (missing(lwid) || is.null(lwid)) 
    #     lwid <- c(keysize, 4)
    # if (missing(lmat) || is.null(lmat)) {
    #     lmat <- rbind(4:3, 2:1)
    #     if (!missing(ColSideColors)) {
    #         if (!is.character(ColSideColors) || length(ColSideColors) != 
    #             nc) 
    #             stop("'ColSideColors' must be a character vector of length ncol(x)")
    #         lmat <- rbind(lmat[1, ] + 1, c(NA, 1), lmat[2, ] + 
    #             1)
    #         lhei <- c(lhei[1], 0.2, lhei[2])
    #     }
    #     if (!missing(RowSideColors)) {
    #         if (!is.character(RowSideColors) || length(RowSideColors) != 
    #             nr) 
    #             stop("'RowSideColors' must be a character vector of length nrow(x)")
    #         lmat <- cbind(lmat[, 1] + 1, c(rep(NA, nrow(lmat) - 
    #             1), 1), lmat[, 2] + 1)
    #         lwid <- c(lwid[1], 0.2, lwid[2])
    #     }
    #     lmat[is.na(lmat)] <- 0
    # }
    # if (length(lhei) != nrow(lmat)) 
    #     stop("lhei must have length = nrow(lmat) = ", nrow(lmat))
    # if (length(lwid) != ncol(lmat)) 
    #     stop("lwid must have length = ncol(lmat) =", ncol(lmat))
    # op <- par(no.readonly = TRUE)
    # on.exit(par(op))
    # layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
    if (!missing(RowSideColors)) {
        par(mar = c(margins[1], 0, 0, 0.5))
        image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    }
    if (!missing(ColSideColors)) {
        par(mar = c(0.5, 0, 0, margins[2]))
        image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    }
    par(mar = c(margins[1], 0, 0, margins[2]))
    x <- t(x)
    cellnote <- t(cellnote)
    if (revC) {
        iy <- nr:1
        if (exists("ddr")) 
            ddr <- rev(ddr)
        x <- x[, iy]
        cellnote <- cellnote[, iy]
    }
    else iy <- 1:nr
    image(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + 
        c(0, nr), axes = FALSE, xlab = "", ylab = "", col = col, 
        breaks = breaks, ...)
    retval$carpet <- x
    if (exists("ddr")) 
        retval$rowDendrogram <- ddr
    if (exists("ddc")) 
        retval$colDendrogram <- ddc
    retval$breaks <- breaks
    retval$col <- col
    if (!invalid(na.color) & any(is.na(x))) {
        mmat <- ifelse(is.na(x), 1, NA)
        image(1:nc, 1:nr, mmat, axes = FALSE, xlab = "", ylab = "", 
            col = na.color, add = TRUE)
    }
    axis(1, 1:nc, labels = labCol, las = 2, line = -0.5, tick = 0, 
        cex.axis = cexCol)
    if (!is.null(xlab)) 
        mtext(xlab, side = 1, line = margins[1] - 1.25)
    axis(4, iy, labels = labRow, las = 2, line = -0.5, tick = 0, 
        cex.axis = cexRow)
    if (!is.null(ylab)) 
        mtext(ylab, side = 4, line = margins[2] - 1.25)
    if (!missing(add.expr)) 
        eval(substitute(add.expr))
    if (!missing(colsep)) 
        for (csep in colsep) rect(xleft = csep + 0.5, ybottom = rep(0, 
            length(csep)), xright = csep + 0.5 + sepwidth[1], 
            ytop = rep(ncol(x) + 1, csep), lty = 1, lwd = 1, 
            col = sepcolor, border = sepcolor)
    if (!missing(rowsep)) 
        for (rsep in rowsep) rect(xleft = 0, ybottom = (ncol(x) + 
            1 - rsep) - 0.5, xright = nrow(x) + 1, ytop = (ncol(x) + 
            1 - rsep) - 0.5 - sepwidth[2], lty = 1, lwd = 1, 
            col = sepcolor, border = sepcolor)
    min.scale <- min(breaks)
    max.scale <- max(breaks)
    x.scaled <- scale01(t(x), min.scale, max.scale)
    if (trace %in% c("both", "column")) {
        retval$vline <- vline
        vline.vals <- scale01(vline, min.scale, max.scale)
        for (i in colInd) {
            if (!is.null(vline)) {
                abline(v = i - 0.5 + vline.vals, col = linecol, 
                  lty = 2)
            }
            xv <- rep(i, nrow(x.scaled)) + x.scaled[, i] - 0.5
            xv <- c(xv[1], xv)
            yv <- 1:length(xv) - 0.5
            lines(x = xv, y = yv, lwd = 1, col = tracecol, type = "s")
        }
    }
    if (trace %in% c("both", "row")) {
        retval$hline <- hline
        hline.vals <- scale01(hline, min.scale, max.scale)
        for (i in rowInd) {
            if (!is.null(hline)) {
                abline(h = i + hline, col = linecol, lty = 2)
            }
            yv <- rep(i, ncol(x.scaled)) + x.scaled[i, ] - 0.5
            yv <- rev(c(yv[1], yv))
            xv <- length(yv):1 - 0.5
            lines(x = xv, y = yv, lwd = 1, col = tracecol, type = "s")
        }
    }
    if (!missing(cellnote)) 
        text(x = c(row(cellnote)), y = c(col(cellnote)), labels = c(cellnote), 
            col = notecol, cex = notecex)
    par(mar = c(margins[1], 0, 0, 0))
    # if (dendrogram %in% c("both", "row")) {
    #     plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
    # }
    # else plot.new()
    par(mar = c(0, 0, if (!is.null(main)) 5 else 0, margins[2]))
    # if (dendrogram %in% c("both", "column")) {
    #     plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
    # }
    # else plot.new()
    if (!is.null(main)) 
        title(main, cex.main = 1.5 * op[["cex.main"]])
    # if (key) {
    #     par(mar = c(5, 4, 2, 1), cex = 0.75)
    #     tmpbreaks <- breaks
    #     if (symkey) {
    #         max.raw <- max(abs(c(x, breaks)), na.rm = TRUE)
    #         min.raw <- -max.raw
    #         tmpbreaks[1] <- -max(abs(x), na.rm = TRUE)
    #         tmpbreaks[length(tmpbreaks)] <- max(abs(x), na.rm = TRUE)
    #     }
    #     else {
    #         min.raw <- min(x, na.rm = TRUE)
    #         max.raw <- max(x, na.rm = TRUE)
    #     }
    #     z <- seq(min.raw, max.raw, length = length(col))
    #     image(z = matrix(z, ncol = 1), col = col, breaks = tmpbreaks, 
    #         xaxt = "n", yaxt = "n")
    #     par(usr = c(0, 1, 0, 1))
    #     lv <- pretty(breaks)
    #     xv <- scale01(as.numeric(lv), min.raw, max.raw)
    #     axis(1, at = xv, labels = lv)
    #     if (scale == "row") 
    #         mtext(side = 1, "Row Z-Score", line = 2)
    #     else if (scale == "column") 
    #         mtext(side = 1, "Column Z-Score", line = 2)
    #     else mtext(side = 1, "Value", line = 2)
    #     if (density.info == "density") {
    #         dens <- density(x, adjust = densadj, na.rm = TRUE)
    #         omit <- dens$x < min(breaks) | dens$x > max(breaks)
    #         dens$x <- dens$x[-omit]
    #         dens$y <- dens$y[-omit]
    #         dens$x <- scale01(dens$x, min.raw, max.raw)
    #         lines(dens$x, dens$y/max(dens$y) * 0.95, col = denscol, 
    #             lwd = 1)
    #         axis(2, at = pretty(dens$y)/max(dens$y) * 0.95, pretty(dens$y))
    #         title("Color Key\nand Density Plot")
    #         par(cex = 0.5)
    #         mtext(side = 2, "Density", line = 2)
    #     }
    #     else if (density.info == "histogram") {
    #         h <- hist(x, plot = FALSE, breaks = breaks)
    #         hx <- scale01(breaks, min.raw, max.raw)
    #         hy <- c(h$counts, h$counts[length(h$counts)])
    #         lines(hx, hy/max(hy) * 0.95, lwd = 1, type = "s", 
    #             col = denscol)
    #         axis(2, at = pretty(hy)/max(hy) * 0.95, pretty(hy))
    #         title("Color Key\nand Histogram")
    #         par(cex = 0.5)
    #         mtext(side = 2, "Count", line = 2)
    #     }
    #     else title("Color Key")
    # }
    # else plot.new()
    retval$colorTable <- data.frame(low = retval$breaks[-length(retval$breaks)], 
        high = retval$breaks[-1], color = retval$col)
    invisible(retval)
}

#rm(list=ls())
param <- list()
param$debug <- FALSE #T

print("Finished loading functions")


if(param$debug) {
    param <- LoadDebugParams(param)
} else {
    param <- LoadReqiredParams(param)
    param <- LoadOptionalParams(param)
}
print (param)
Rprof() #Remove me
# Load Chip data, alread in pval format
# Column.order need to be a vecor of strings, that correspond to the order (and
# inclustion) of chip experiments.
chip <- GetChipData(param$annotated.macs.file, column.order = param$chip.order)

#Get the top percentages on different criteria
#ix.top <- GetTopRowsFromMatrix(chip$peaks, percentage = param$filter.percentage)
#chip$peaks <- chip$peaks[ix.top$union, ]
#chip$targets <- chip$targets[ix.top$union]
#Bin data:
chip$peaks <- bindata.non.zero.matrix(chip$peaks, qunts = seq(0, 0.9, length=(param$number.bins+1)))

file.names <- unlist(strsplit(param$rna.files, split="::"))
print(file.names)
file.lables <- unlist(strsplit(param$rna.names, split="::"))
print(file.lables)

library("RColorBrewer")
color.2 =  rev(colorRampPalette(c("darkblue", "steelblue", "lightgrey"))(param$number.bins))
if(!is.na(file.names) && file.names != "none") {
    #rna.scores <- GetRNAData(file.names, file.lables = file.lables, fpkm = "avg")
    #do differential expression if the user specifies such:
    if(param$rna.normalization != "no") {
        if(param$rna.normalization != "mean"){
            norm.file.names <- unlist(strsplit(param$rna.normalization.file, split="::"))
            all.file.names <- c(file.names, norm.file.names)
            all.file.lables <- c(file.lables, norm.file.names)
            print (all.file.names)
            print (all.file.lables)
            rna.scores <- GetRNAData(all.file.names, file.lables = all.file.lables, fpkm = "avg")
            rna.scores <- NormalizeRNA(scores = rna.scores)
            rna.scores.sign <- rna.scores/abs(rna.scores)
            rna.scores.sign[which(is.nan(rna.scores.sign))] <- 0
            rna.scores <- bindata.non.zero.matrix(abs(rna.scores), qunts = seq(0, 0.9, length=(param$number.bins+1)))
            rna.scores <- rna.scores.sign * rna.scores
            color.2 = colorRampPalette(c("darkblue", "steelblue", "lightgrey", "pink", "darkred"))(param$number.bins)
        } else {
            print("we are normalizing by average")
            print("file lables are:")
            print(file.lables)
            rna.scores <- GetRNAData(file.names, file.lables = file.lables, fpkm = "avg")
            rna.scores <- t(apply(rna.scores,1, function(x) { if(mean(x)!=0){ return(x/mean(x)) } else { return(x) }  }))
            rna.scores <- bindata.non.zero.matrix(rna.scores, qunts = seq(0, 0.9, length=(param$number.bins+1)))
        }
    } else {
        rna.scores <- GetRNAData(file.names, file.lables = file.lables, fpkm = "avg")
        rna.scores <- bindata.non.zero.matrix(rna.scores, qunts = seq(0, 0.9, length=(param$number.bins+1)))
    }
    
    # chip$peaks <- chip$peaks[order(chip$targets), ]
    # rna.scores <- rna.scores[order(rownames(rna.scores)), ]
    # chip$targets <- chip$targets[order(chip$targets)]

    rna.scores.mapped <- MapExpressiontoChip(chip.targets = chip$targets, expression=rna.scores)
    all.data <- cbind(chip$peaks, rna.scores.mapped)
    if(param$include.targetless!="yes"){
        all.data <- all.data[-which(chip$targets==""),]
        chip$peaks <- chip$peaks[-which(chip$targets==""),]
        rna.scores.mapped <- rna.scores.mapped[-which(chip$targets==""),]
        chip$targets <- chip$targets[-which(chip$targets=="")]
    }
} else {
    all.data <- chip$peaks
}

print("Clustering")
set.seed(1234)
km <- kmeans(all.data, param$clustering.number.of.clusters, iter.max = 50)
print("Ordering")
km.order <- GenerateKMOrder(all.data, km)

##AM edits##
km.new.order <- numeric(length(km$cluster))
for(i in 1:length(km.order)){
    cur.clst <- km.order[i]
    ix <- which(km$cluster==cur.clst)
    km.new.order[ix] <- i
}
km.order <- 1:length(km.order)
km$cluster <- km.new.order
## Done AM ##
print("Creating image")
#bmp(param$heatmap.document.name, 640, 480)#1280, 960)
bmp(paste(param$heatmap.document.name,".bmp", sep = ""), 5120, 3840)#2560, 1920)
#pdf(paste(param$heatmap.document.name,".pdf", sep = ""))
color.1 <- c("#000000", colorRampPalette(c("blue", "yellow","orange","red"))(param$number.bins))
if(!is.na(file.names) && file.names != "none") {
    PrintClusters(trgts=chip$targets,
        k.ix=km$cluster,
        f.nm = param$cluster.groups.document.name,
        km.order = NA)

    #if(param$rna.normalization != "no") {
    #    data.split <- cbind(chip$peaks, PrepareRNAforHeatmap(rna.scores.mapped))
    #} else {
    #    data.split <- cbind(chip$peaks, rna.scores.mapped)
    #}
    expression.width.multiplier <- 2
    if(ncol(rna.scores.mapped)==1) {
        rna.scores.mapped <- cbind(rna.scores.mapped, rna.scores.mapped)
        expression.width.multiplier <- 1
    }

    layout(matrix(c(1,2,2,5,
        1,3,4,6,
        1,3,4,7,
        1,3,4,8,
        1,3,4,9),
    5,4,
    byrow=T),
    widths=c(1,2*ncol(chip$peaks),expression.width.multiplier*ncol(rna.scores.mapped),1),
    heights=c(1,10,2,10,2))
    #1
    plot.new()
    #2
    plot.new()
    #3
    print("Creating peak image")
    CreateIndividualHeatMap(chip$peaks,
        km,
        km.order, color.ramp=color.1)
    #4
    print("Creating rna image")
    CreateIndividualHeatMap(rna.scores.mapped,
        km,
        km.order, color.ramp=color.2)
    #5
    plot.new()
    #6
    image(t(matrix(1:param$number.bins,nc=1)), col=color.1)
    #7
    plot.new()
    #8
    image(t(matrix(1:param$number.bins,nc=1)), col=color.2)
    #9
    plot.new()
} else {
    PrintClusters(trgts=1:nrow(chip$peaks),
        k.ix=km$cluster,
        f.nm = param$cluster.groups.document.name,
        km.order = NA)


    layout(matrix(c(1,2,2,
        1,3,4,
        1,3,5),3,3,byrow=T),
        widths=c(1,3*ncol(chip$peaks),1),
        heights=c(1,8,1))
    #1
    plot.new()
    #2
    plot.new()
    #3
    CreateIndividualHeatMap(chip$peaks,
        km,
        km.order, color.ramp=color.1)
    #4
    image(t(matrix(1:param$number.bins,nc=1)), col=color.1)
    #5
    plot.new()

}

dev.off()
Rprof(NULL)
profile <- summaryRprof()
print(str(profile))
print(profile)
#save.image("test/testData.RData")