--- a/create_heatmap.R	Fri Dec 15 16:48:53 2017 -0500
+++ b/create_heatmap.R	Mon Dec 18 14:30:14 2017 -0500
@@ -3,8 +3,8 @@
 suppressPackageStartupMessages(library("optparse"))
 
 option_list <- list(
-            make_option(c("-i", "--input_dir"), action="store", dest="input_dir", help="IDEAS para files directory"),
-            make_option(c("-o", "--output_dir"), action="store", dest="output_dir", help="PDF output directory")
+        make_option(c("-i", "--input_dir"), action="store", dest="input_dir", help="IDEAS para files directory"),
+        make_option(c("-o", "--output_dir"), action="store", dest="output_dir", help="PDF output directory")
 )
 
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
@@ -16,98 +16,141 @@
     # based on the received data_frame which was created
     # by reading the .para file.
     num_columns = dim(data_frame)[2];
+    cat("num_columns: ", num_columns, "\n");
     num_rows = dim(data_frame)[1];
-    p = (sqrt(9 + 8 * (num_columns - 1)) - 3) / 2;
+    cat("num_rows: ", num_rows, "\n");
+    p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
+    cat("9 + 8 * (num_columns-1): ", 9 + 8 * (num_columns-1), "\n");
+    cat("sqrt(9 + 8 * (num_columns-1)): ", sqrt(9 + 8 * (num_columns-1)), "\n");
+    cat("p: ", p, "\n");
     data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
+    cat("dim(data_matrix)[1]: ", dim(data_matrix)[1], "\n");
+    cat("dim(data_matrix)[2]: ", dim(data_matrix)[2], "\n");
     colnames(data_matrix) = colnames(data_frame)[1+1:p];
-    marks = colnames(data_matrix);
-    rownames(data_matrix) = paste(1:num_rows, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
-
+    cat("colnames(data_matrix): ", colnames(data_matrix), "\n");
+    histone_marks = colnames(data_matrix);
+    cat("histone_marks: ", histone_marks, "\n");
+    rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
+    cat("rownames(data_matrix): ", rownames(data_matrix), "\n");
     # Open the output PDF file.
     pdf(file=output_file_name);
     # Set graphical parameters.
     par(mar=c(6, 1, 1, 6));
     # Create a vector containing the minimum and maximum values in data_matrix.
     min_max_vector = range(data_matrix);
-    # Define colors for the palette.
-    colors = 0:100 / 100 * (min_max_vector[2] - min_max_vector[1]) + min_max_vector[1];
-    # Create the color palette.
+    cat("min_max_vector: ", min_max_vector, "\n");
+    # Create a color palette.
     my_palette = colorRampPalette(c("white", "dark blue"))(n=100);
     defpalette = palette(my_palette);
     # Plot the heatmap for the current .para / .state combination.
     plot(NA, NA, xlim=c(0, p+0.7), ylim=c(0, num_rows), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
     axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
     axis(4, at=1:num_rows-0.5, labels=rownames(data_matrix), las=2);
-    rect(rep(1:p-1, num_rows), rep(1:num_rows, each=p), rep(1:p, num_rows), rep(1:num_rows, each=p),
-            col=round((t(data_matrix)-min_max_vector[1])/(min_max_vector[2]-min_max_vector[1])*100));
-    markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
+    color = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100);
+    cat("color: ", color, "\n");
+    rect(rep(1:p-1, num_rows), rep(1:num_rows-1, each=p), rep(1:p, num_rows), rep(1:num_rows, each=p), col=color);
+    histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
+    cat("histone_mark_color: ", histone_mark_color, "\n");
 
-    for(i in 1:length(marks)) {
-        if(regexpr("h3k4me3", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(255, 0, 0);
-        }
-        if(regexpr("h3k4me2", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(250, 100, 0);
+    # Specify a color for common feature names like "h3k4me3".
+    # These are histone marks frequently used to identify
+    # promoter activities in a cell, and are often displayed
+    # in shades of red.
+    for(i in 1:length(histone_marks)) {
+        cat("histone_marks[i]: ", histone_marks[i], "\n");
+        if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(255, 0, 0);
+            cat("histone_mark_color[i]: ", histone_mark_color[i], "\n");
         }
-        if(regexpr("h3k4me1", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(250, 250, 0);
+        if(regexpr("h3k4me2", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 100, 0);
         }
-        if(regexpr("h3k36me3", tolower(marks[i]))>0) {
-            markcolor[i,] = c(0, 150, 0);
+        if(regexpr("h3k4me1", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 250, 0);
         }
-        if(regexpr("h2a", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(0, 150, 150);
+        if(regexpr("h3k36me3", tolower(histone_marks[i]))>0) {
+            histone_mark_color[i,] = c(0, 150, 0);
         }
-        if(regexpr("dnase", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(0, 200, 200);
+        if(regexpr("h2a", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 150, 150);
         }
-        if(regexpr("h3k9ac", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(250, 0, 200);
+        if(regexpr("dnase", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 200, 200);
         }
-        if(regexpr("h3k9me3", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(100, 100, 100);
+        if(regexpr("h3k9ac", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 0, 200);
         }
-        if(regexpr("h3k27ac", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(250, 150, 0);
+        if(regexpr("h3k9me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(100, 100, 100);
         }
-        if(regexpr("h3k27me3", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(0, 0, 200);
+        if(regexpr("h3k27ac", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 150, 0);
         }
-        if(regexpr("h3k79me2", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(200, 0, 200);
+        if(regexpr("h3k27me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 0, 200);
         }
-        if(regexpr("h4k20me1", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(50, 200, 50);
+        if(regexpr("h3k79me2", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(200, 0, 200);
         }
-        if(regexpr("ctcf", tolower(marks[i])) > 0) {
-            markcolor[i,] = c(200, 0, 250);
+        if(regexpr("h4k20me1", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(50, 200, 50);
         }
-        state_color = get_state_color(data_matrix, markcolor)[,2];
+        if(regexpr("ctcf", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(200, 0, 250);
+        }
+        state_color = get_state_color(data_matrix, histone_mark_color)[,2];
+        cat("state_color: ", state_color, "\n");
     }
-    rect(rep(p+0.2,num_rows), 1:num_rows-0.8, rep(p+0.8,num_rows), 1:num_rows-0.2, col=state_color);
+    rect(rep(p+0.2, num_rows), 1:num_rows-0.8, rep(p+0.8, num_rows), 1:num_rows-0.2, col=state_color);
     palette(defpalette);
     dev.off();
 }
 
-get_state_color<-function(statemean, markcolor) {
-    rg=apply(statemean, 1, range);
-    mm=NULL;
-    for(i in 1:dim(statemean)[1]) {
-        mm = rbind(mm, (statemean[i,]-rg[1,i]+1e-10)/(rg[2,i]-rg[1,i]+1e-10));
+get_state_color <- function(data_matrix, histone_mark_color) {
+    cat("XXX histone_mark_color: ", histone_mark_color, "\n");
+    #histone_mark_color = rep("", dim(data_matrix)[2]);
+    #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
+    #cat("XXX order(apply(data_matrix, 2, sd), decreasing=T): ", order(apply(data_matrix, 2, sd), decreasing=T), "\n");
+    #cat("XXX dim(data_matrix)[2]-1: ", dim(data_matrix)[2]-1, "\n");
+    #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
+    #histone_mark_color[order(apply(data_matrix, 2, sd), decreasing=T)] = hsv((1:dim(data_matrix)[2]-1) / dim(data_matrix)[2], 1, 1);
+    #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
+    #histone_mark_color = t(col2rgb(histone_mark_color));
+    #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
+    rg = apply(data_matrix, 1, range);
+    cat("XXX rg: ", rg, "\n");
+    mm = NULL;
+    cat("XXX dim(data_matrix): ", dim(data_matrix), "\n");
+    cat("XXX dim(data_matrix)[1]: ", dim(data_matrix)[1], "\n");
+    for(i in 1:dim(data_matrix)[1]) {
+        mm = rbind(mm, (data_matrix[i,] - rg[1, i] + 1e-10) / (rg[2, i] -rg[1, i] + 1e-10));
+        cat("XXX mm: ", mm, "\n");
     }
     mm = mm^5;
+    cat("XXX dim(mm): ", dim(mm), "\n");
+    cat("XXX dim(mm)[2]: ", dim(mm)[2], "\n");
     if(dim(mm)[2] > 1) {
         mm = mm / (apply(mm, 1, sum) + 1e-10);
+        cat("XXX mm: ", mm, "\n");
     }
-    mycol = mm%*%markcolor;
-    s = apply(statemean, 1, max);
-    s = (s-min(s)) / (max(s)-min(s) + 1e-10);
-    mycol = round(255-(255-mycol)*s/0.5);
-    mycol[mycol<0] = 0;
-    rt = paste(mycol[,1], mycol[,2], mycol[,3], sep=",");
-    h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
+    state_color = mm%*%histone_mark_color;
+    cat("XXX state_color: ", state_color, "\n");
+    s = apply(data_matrix, 1, max);
+    cat("XXX s: ", s, "\n");
+    s = (s - min(s)) / (max(s) - min(s) + 1e-10);
+    cat("XXX s: ", s, "\n");
+    state_color = round(255 - (255 - state_color) * s/0.5);
+    cat("XXX state_color: ", state_color, "\n");
+    state_color[state_color<0] = 0;
+    cat("XXX state_color: ", state_color, "\n");
+    rt = paste(state_color[,1], state_color[,2], state_color[,3], sep=",");
+    cat("XXX rt: ", rt, "\n");
+    h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
+    cat("XXX h: ", h, "\n");
     h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
+    cat("XXX h: ", h, "\n");
     rt = cbind(rt, h);
+    cat("XXX rt: ", rt, "\n");
     return(rt);
 }