Mercurial > repos > greg > ideas_genome_tracks

changeset 47:c31d53123fae draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 14 Dec 2017 13:35:49 -0500
parents 1290764e2750
children 559971a4f55e
files ideas_genome_tracks.R
diffstat 1 files changed, 26 insertions(+), 119 deletions(-) [+]
line wrap: on
line diff
--- a/ideas_genome_tracks.R	Thu Dec 14 13:35:41 2017 -0500
+++ b/ideas_genome_tracks.R	Thu Dec 14 13:35:49 2017 -0500
@@ -2,7 +2,7 @@
 
 suppressPackageStartupMessages(library("data.table"))
 suppressPackageStartupMessages(library("optparse"))
-suppressPackageStartupMessages(library("viridislite"))
+suppressPackageStartupMessages(library("viridisLite"))
 
 option_list <- list(
         make_option(c("--build"), action="store", dest="build", help="Genome build"),
@@ -10,7 +10,6 @@
         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"),
@@ -35,9 +34,9 @@
     cat(s, file=output_trackhub);
 }
 
-create_track = function(input_dir_para, input_dir_state, chrom_len_file, base_track_file_name) {
-    para_files <- list.files(path=input_dir_para, full.names=TRUE);
-    cat("para_files: ", para_files, "\n");
+create_track = function(input_dir_state, chrom_len_file, base_track_file_name) {
+    # Create everythin 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);
     g = NULL;
@@ -71,14 +70,11 @@
     L = dim(g1)[1];
     n = dim(state)[2];
     cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
-    g1 = NULL;
     options(scipen=999);
     tt = which(chr[2:L]!=chr[2:L-1]);
     tt = c(tt, which(posst[2:L]!=posed[2:L-1]));
     tt = sort(unique(tt));
-    state_color <- get_state_color(para_files[1])
     for(i in 1:n) {
-        cat("i: ", i, "\n");
         tstate = state[,i];
         t = c(tt, which(tstate[2:L]!=tstate[2:L-1]));
         t = sort(unique(t));
@@ -87,7 +83,7 @@
         np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
         track_file_name_bed <- get_track_file_name(base_track_file_name, i, "bed");
         track_file_name_bigbed <- get_track_file_name(base_track_file_name, i, "bigbed");
-        x = cbind(np[, 1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]);
+        x = cbind(np[, 1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3]);
         write.table(as.matrix(x), track_file_name_bed, quote=F, row.names=F, col.names=F);
         system(paste("bedToBigBed ", track_file_name_bed, chrom_len_file, " ", track_file_name_bigbed));
         system(paste("rm ", track_file_name_bed));
@@ -95,30 +91,37 @@
     return(cells);
 }
 
-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) {
+create_track_db = function(galaxy_url, encoded_dataset_id, input_dir_state, build, chrom_len_file, tracks_dir, hub_name, short_label, long_label) {
+    # Create a trackDb.txt file for each state.
     base_track_file_name <- paste(tracks_dir, hub_name, sep="");
-    cells = create_track(input_dir_para, input_dir_state, chrom_len_file, base_track_file_name);
-    cell_info = cbind(cells, cells, cells, "#000000");
-    cell_info = array(cell_info, dim=c(length(cells), 4));
-    cell_info = as.matrix(cell_info);
+    cells = create_track(input_dir_state, chrom_len_file, base_track_file_name);
+    # Create a a character vector of 1024 viridis color hex codes.
+    viridis_vector <- viridis(1024, alpha=1, begin=0, end=1, direction=1, option="D");
+    colors_used <- vector();
     track_db = NULL;
     for (i in 1:length(cells)) {
+        # 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.
+        viridis_color_hex_code <- viridis_vector[color_index]
         big_data_url <- get_big_data_url(galaxy_url, encoded_dataset_id, tracks_dir, i, build);
-        ii = which(cells[i] == cell_info[,1]);
-        if (length(ii) == 0) {
-            next;
-        }
-        ii = ii[1];
         # trackDb.txt track entry.
         track_db = c(track_db, paste("hub ", hub_name, "_track_", i, sep=""));
         track_db = c(track_db, "type bigBed");
         track_db = c(track_db, paste("bigDataUrl", big_data_url, sep=" "));
         track_db = c(track_db, paste("shortLabel", short_label, sep=" "));
         track_db = c(track_db, paste("longLabel", long_label, sep=" "));
-        track_db = c(track_db, paste("priority", ii));
+        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(cell_info[ii,4])), collapse=","), sep=" "));
+        track_db = c(track_db, paste("color", paste(c(col2rgb(viridis_color_hex_code)), collapse=","), sep=" "));
         track_db = c(track_db, "visibility dense");
         track_db = c(track_db, "");
     }
@@ -131,102 +134,6 @@
     return(s)
 }
 
-get_rgb<-function(statemean, markcolor=NULL)
-{
-    if(length(markcolor) == 0) {
-        markcolor = rep("",dim(statemean)[2]);
-        markcolor[order(apply(statemean,2,sd),decreasing=T)]=hsv((1:dim(statemean)[2]-1)/dim(statemean)[2],1,1)
-        markcolor = t(col2rgb(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));
-    }
-    mm = mm^5;
-    if (dim(mm)[2]>1) {
-        mm = mm / (apply(mm, 1, sum)+1e-10);
-    }
-    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])}));
-    h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
-    rt = cbind(rt, h);
-    return(rt);
-
-    h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
-    h[,2] = h[,2]*s;
-    h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
-    rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}), h);
-    return(rt);
-}
-
-get_state_color <- function(para_file, cols=c("white", "dark blue")) {
-    x = read.table(para_file, comment="!", header=T);
-    k = dim(x)[2];
-    l = dim(x)[1];
-    p = (sqrt(9+8*(k-1))-3)/2;
-    m = as.matrix(x[,1+1:p]/x[,1]);
-    colnames(m) = colnames(x)[1+1:p];
-    marks = colnames(m);
-    rownames(m) = paste(1:l-1," (",round(x[,1]/sum(x[,1])*10000)/100,"%)",sep="");
-    par(mar=c(6,1,1,6));
-    rg = range(m);
-    colors = 0:100/100*(rg[2]-rg[1])+rg[1];
-    markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
-    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);
-        }
-        if (regexpr("h3k4me1",tolower(marks[i]))>0) {
-            markcolor[i,]=c(250,250,0);
-        }
-        if (regexpr("h3k36me3",tolower(marks[i]))>0) {
-            markcolor[i,]=c(0,150,0);
-        }
-        if (regexpr("h2a",tolower(marks[i]))>0) {
-            markcolor[i,]=c(0,150,150);
-        }
-        if (regexpr("dnase",tolower(marks[i]))>0) {
-            markcolor[i,]=c(0,200,200);
-        }
-        if (regexpr("atac",tolower(marks[i]))>0) {
-            markcolor[i,]=c(0,200,200);
-        }
-        if (regexpr("h3k9ac",tolower(marks[i]))>0) {
-            markcolor[i,]=c(250,0,200);
-        }
-        if (regexpr("h3k9me3",tolower(marks[i]))>0) {
-            markcolor[i,]=c(100,100,100);
-        }
-        if (regexpr("h3k27ac",tolower(marks[i]))>0) {
-            markcolor[i,]=c(250,150,0);
-        }
-        if (regexpr("h3k27me3",tolower(marks[i]))>0) {
-            markcolor[i,]=c(0,0,225);
-        }
-        if (regexpr("h3k79me2",tolower(marks[i]))>0) {
-            markcolor[i,]=c(200,0,200);
-        }
-        if (regexpr("h4k20me1",tolower(marks[i]))>0) {
-            markcolor[i,]=c(50,200,50);
-        }
-        if (regexpr("ctcf",tolower(marks[i]))>0) {
-            markcolor[i,]=c(200,0,250);
-        }
-    }
-    statecolor = get_rgb(m, markcolor)[,];
-    return(statecolor);
-}
-
 get_track_file_name = function(base_track_file_name, index, ext) {
     track_file_name <- paste(base_track_file_name, index, ext, sep=".");
     return(track_file_name);
@@ -254,11 +161,11 @@
 # Create the tracks.
 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_para, opt$input_dir_state, opt$build, opt$chrom_len_file, tracks_dir, opt$hub_name, opt$short_label, opt$long_label);
+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);
 
 # Create the trackDb.txt output.
 track_db_file_path <- paste(tracks_dir, "trackDb.txt", sep="");
 write.table(track_db, file=track_db_file_path, quote=F, row.names=F, col.names=F);
 
 # Create the primary HTML dataset.
-create_primary_html(opt$output_trackhub, tracks_dir, opt$build);
\ No newline at end of file
+create_primary_html(opt$output_trackhub, tracks_dir, opt$build);