# HG changeset patch
# User greg
# Date 1513787164 18000
# Node ID f0c488b73dbc65da22a0768aeee28e1daffe819b
# Parent 794d2104f83d631affa26b69a993fa3f75647672
Uploaded
diff -r 794d2104f83d -r f0c488b73dbc create_heatmap.R
--- /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);
+}
diff -r 794d2104f83d -r f0c488b73dbc ideas_genome_tracks.R
--- 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="");
diff -r 794d2104f83d -r f0c488b73dbc ideas_genome_tracks.xml
--- 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 @@
bedops
r-data.table
r-optparse
- r-viridislite
ucsc-bedtobigbed
-
+