Mercurial > repos > greg > ideas_genome_tracks

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/create_heatmap.R	Wed Dec 20 11:26:04 2017 -0500
@@ -0,0 +1,109 @@
+#!/usr/bin/env Rscript
+
+create_heatmap<-function(data_frame, output_file_name=NULL) {
+    # Plot a heatmap for a .para / .state combination
+    # based on the received data_frame which was created
+    # by reading the .para file.
+    num_columns = dim(data_frame)[2];
+    num_rows = dim(data_frame)[1];
+    p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
+    data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
+    colnames(data_matrix) = colnames(data_frame)[1+1:p];
+    histone_marks = colnames(data_matrix);
+    rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
+    if (!is.null(output_file_name)) {
+        # 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);
+    # 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);
+    color = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100);
+    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]));
+
+    # 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)) {
+        if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(255, 0, 0);
+        }
+        if(regexpr("h3k4me2", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 100, 0);
+        }
+        if(regexpr("h3k4me1", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 250, 0);
+        }
+        if(regexpr("h3k36me3", tolower(histone_marks[i]))>0) {
+            histone_mark_color[i,] = c(0, 150, 0);
+        }
+        if(regexpr("h2a", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 150, 150);
+        }
+        if(regexpr("dnase", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 200, 200);
+        }
+        if(regexpr("h3k9ac", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 0, 200);
+        }
+        if(regexpr("h3k9me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(100, 100, 100);
+        }
+        if(regexpr("h3k27ac", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(250, 150, 0);
+        }
+        if(regexpr("h3k27me3", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(0, 0, 200);
+        }
+        if(regexpr("h3k79me2", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(200, 0, 200);
+        }
+        if(regexpr("h4k20me1", tolower(histone_marks[i])) > 0) {
+            histone_mark_color[i,] = c(50, 200, 50);
+        }
+        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];
+    }
+    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);
+    if (!is.null(output_file_name)) {
+        dev.off();
+    }
+    return(state_color);
+}
+
+get_state_color <- function(data_matrix, histone_mark_color) {
+    range_vector = apply(data_matrix, 1, range);
+    mm = NULL;
+    for(i in 1:dim(data_matrix)[1]) {
+        range_val1 = range_vector[1, i] + 1e-10;
+        range_val2 = range_vector[2, i];
+        mm = rbind(mm, (data_matrix[i,] - range_val1) / (range_val2 - range_val1));
+    }
+    mm = mm^5;
+    if(dim(mm)[2] > 1) {
+        mm = mm / (apply(mm, 1, sum) + 1e-10);
+    }
+    state_color = mm%*%histone_mark_color;
+    s = apply(data_matrix, 1, max);
+    s = (s - min(s)) / (max(s) - min(s) + 1e-10);
+    state_color = round(255 - (255 - state_color) * s/0.5);
+    state_color[state_color<0] = 0;
+    rt = paste(state_color[,1], state_color[,2], state_color[,3], sep=",");
+    h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
+    h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
+    rt = cbind(rt, h);
+    return(rt);
+}
--- a/ideas_genome_tracks.R	Tue Dec 19 12:53:24 2017 -0500
+++ b/ideas_genome_tracks.R	Wed Dec 20 11:26:04 2017 -0500
@@ -2,7 +2,6 @@

 suppressPackageStartupMessages(library("data.table"))
 suppressPackageStartupMessages(library("optparse"))
-suppressPackageStartupMessages(library("viridisLite"))

 option_list <- list(
     make_option(c("--build"), action="store", dest="build", help="Genome build"),
@@ -10,6 +9,7 @@
     make_option(c("--email"),  action="store", dest="email", help="User email address"),
     make_option(c("--galaxy_url"),  action="store", dest="galaxy_url", help="Galaxy instance base URL"),
     make_option(c("--hub_name"),  action="store", dest="hub_name", default=NULL, help="Hub name without spaces"),
+    make_option(c("--input_dir_para"), action="store", dest="input_dir_para", help="Directory containing .para outputs from IDEAS"),
     make_option(c("--input_dir_state"), action="store", dest="input_dir_state", help="Directory containing .state outputs from IDEAS"),
     make_option(c("--long_label"), action="store", dest="long_label", help="Hub long label"),
     make_option(c("--output_trackhub"),  action="store", dest="output_trackhub", help="Output hub file"),
@@ -37,7 +37,7 @@
 }

 create_track = function(input_dir_state, chrom_len_file, base_track_file_name) {
-    # Create everythin needed, including the bigbed file,
+    # Create everything needed, including the bigbed file,
     # to render the tracks within the UCSC track hub.
     state_files <- list.files(path=input_dir_state, full.names=TRUE);
     genome_size = read.table(chrom_len_file);
@@ -98,15 +98,12 @@
     return(cells);
 }

-create_track_db = function(galaxy_url, encoded_dataset_id, input_dir_state, build, chrom_len_file, tracks_dir, hub_name, short_label, long_label, state_indexes, state_colors) {
+create_track_db = function(galaxy_url, encoded_dataset_id, input_dir_para, input_dir_state, build,
+        chrom_len_file, tracks_dir, hub_name, short_label, long_label, state_indexes, state_colors) {
     # Create a trackDb.txt file that includes each state.
+    para_files <- list.files(path=input_dir_para, full.names=TRUE);
     base_track_file_name <- paste(tracks_dir, hub_name, sep="");
     cells = create_track(input_dir_state, chrom_len_file, base_track_file_name);
-    # Create a a character vector of 1024 viridis color hex codes.
-    # This vector could be used even if the state_colors were received
-    # since colors may not have been chosen for all states.
-    viridis_vector <- viridis(1024, alpha=1, begin=0, end=1, direction=1, option="D");
-    colors_used <- vector();
     if (!is.null(state_indexes)) {
         # Split state_indexes into a list of integers.
         s_indexes <- c();
@@ -126,20 +123,12 @@
     track_db = NULL;
     for (i in 1:length(cells)) {
         if (is.null(state_indexes) || !is.element(i, s_indexes)) {
-            # Generate a random number between 1 and 1024 as
-            # the viridis_vector index for the next state color.
-            color_index <- sample(1:1024, 1);
-            # Make sure the color was not previously chosen.
-            while(is.element(color_index, colors_used)) {
-                color_index <- sample(1:1024, 1);
-            }
-            # Append the color to our list of chosen colors.
-            append(colors_used, color_index);
-            # Get the hex code from viridis_vector.
-            color_hex_code <- viridis_vector[color_index];
+            data_frame <- read.table(para_files[i], comment="!", header=T);
+            color <- create_heatmap(data_frame);
         } else {
             # Use the selected color for the current state.
             color_hex_code <- s_colors[i];
+            color <- paste(c(col2rgb(color_hex_code)), collapse=",");
         }
         big_data_url <- get_big_data_url(galaxy_url, encoded_dataset_id, tracks_dir, i, build);
         # trackDb.txt track entry.
@@ -151,7 +140,7 @@
         track_db = c(track_db, paste("priority", i));
         track_db = c(track_db, "itemRgb on");
         track_db = c(track_db, "maxItems 100000");
-        track_db = c(track_db, paste("color", paste(c(col2rgb(color_hex_code)), collapse=","), sep=" "));
+        track_db = c(track_db, paste("color", color, sep=" "));
         track_db = c(track_db, "visibility dense");
         track_db = c(track_db, "");
     }
@@ -189,9 +178,11 @@
 write.table(contents, file=genomes_file_path, quote=F, row.names=F, col.names=F);

 # Create the tracks.
+source("create_heatmap.R");
 tracks_dir <- paste(hub_dir, opt$build, "/", sep="");
 dir.create(tracks_dir, showWarnings=FALSE);
-track_db <- create_track_db(opt$galaxy_url, opt$output_trackhub_id, opt$input_dir_state, opt$build, opt$chrom_len_file, tracks_dir, opt$hub_name, opt$short_label, opt$long_label, opt$state_indexes, opt$state_colors);
+track_db <- create_track_db(opt$galaxy_url, opt$output_trackhub_id, opt$input_dir_state, opt$input_dir_state, opt$build,
+        opt$chrom_len_file, tracks_dir, opt$hub_name, opt$short_label, opt$long_label, opt$state_indexes, opt$state_colors);

 # Create the trackDb.txt output.
 track_db_file_path <- paste(tracks_dir, "trackDb.txt", sep="");
--- a/ideas_genome_tracks.xml	Tue Dec 19 12:53:24 2017 -0500
+++ b/ideas_genome_tracks.xml	Wed Dec 20 11:26:04 2017 -0500
@@ -4,7 +4,6 @@
         <requirement type="package" version="2.4.27">bedops</requirement>
         <requirement type="package" version="1.10.4">r-data.table</requirement>
          <requirement type="package" version="1.4.4">r-optparse</requirement>
-         <requirement type="package" version="0.2.0">r-viridislite</requirement>
          <requirement type="package" version="332">ucsc-bedtobigbed</requirement>
     </requirements>
     <command detect_errors="exit_code"><![CDATA[
@@ -19,6 +18,7 @@
 #else:
     #set email = $__user__.email
 #end if
+#set input_dir_state = "input_dir_para"
 #set input_dir_state = "input_dir_state"
 #set select_state_color = $select_state_color_cond.select_state_color
 #if str($select_state_color) == 'yes':
@@ -32,6 +32,7 @@
     #set state_indexes = ",".join($state_indexes)
     #set state_colors = ",".join($state_colors)
 #end if
+mkdir $input_dir_para &&
 mkdir $input_dir_state &&
 mkdir '$output_trackhub.files_path' &&
 #for $i in $input:
@@ -40,7 +41,9 @@
     #end if
     #set filename = $i.file_name
     #set name = $i.name
-    #if $name.endswith(".state"):
+    #if $name.endswith(".para"):
+        ln -s $filename $input_dir_para/$name &&
+    #else if $name.endswith(".state"):
         ln -s $filename $input_dir_state/$name &&
     #end if
 #end for
@@ -50,6 +53,7 @@
 --email '$email'
 --galaxy_url $galaxy_url
 --hub_name '$hub_name'
+--input_dir_para '$input_dir_para'
 --input_dir_state '$input_dir_state'
 --long_label '$sanitized_long_label'
 --output_trackhub '$output_trackhub'
@@ -76,7 +80,7 @@
             </param>
             <when value="no" />
             <when value="yes">
-                <repeat name="state_color_repeat" title="State colors" min="1">
+                <repeat name="state_color_repeat" title="state color" min="1">
                     <param name="state_index" type="integer" value="" min="1" label="State index" help="Enter the integer index of the desired state">
                         <validator type="empty_field"/>
                     </param>