changeset 5:d1703a15c393 draft

Uploaded

--- a/ideas_genome_tracks.R	Wed Nov 15 13:13:52 2017 -0500
+++ b/ideas_genome_tracks.R	Thu Nov 16 09:04:48 2017 -0500
@@ -1,5 +1,6 @@
 #!/usr/bin/env Rscript
 
+suppressPackageStartupMessages(library("data.table"))
 suppressPackageStartupMessages(library("optparse"))
 
 option_list <- list(
@@ -23,20 +24,21 @@
 
 create_color_scheme = function(input_dir_para) {
     # Create the color scheme.
-    para_files <- list.files(path=input_dir_para, full.names=True);
+    para_files <- list.files(path=input_dir_para, full.names=TRUE);
     mc = NULL;
-    x = read.table(opt$parameter_files, comment="!", nrows=1);
+    x = read.table(para_files[1], comment="!", nrows=1);
     l = as.integer(regexpr("\\*", as.matrix(x)));
     l = min(which(l>0)) - 2;
-    x = as.matrix(read.table(opt$parameter_files));
+    x = as.matrix(read.table(para_files[1]));
     if (length(para_files) > 1) {
         for (i in 2:length(para_files)) {
-            x = x + as.matrix(read.table(paste(para_files[i], ".para", sep="")));
+            x = x + as.matrix(read.table(para_files[i]));
         }
     }
     p = x[,1] / sum(x[,1]);
     m = array(as.matrix(x[,1:l+1] / x[,1]), dim=c(dim(x)[1], l));
-    return get_state_color(m, mc);
+    state_color <- get_state_color(m, mc);
+    state_color
 }
 
 get_state_color = function(statemean, markcolor=NULL)
@@ -56,19 +58,18 @@
     mycol = mm%*%markcolor;
     s = apply(statemean, 1, max);
     s = (s - min(s)) / (max(s) - min(s) + 1e-10);
-    h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
+    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);
+    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);
 }
 
 create_track = function(input_dir_state, chrom_len_file, outpref, state_color, header, state_name) {
-    state_files <- list.files(path=input_dir_state, full.names=True);
+    state_files <- list.files(path=input_dir_state, full.names=TRUE);
     genome_size = read.table(chrom_len_file);
     g = NULL;
     for(i in state_files) {
-        message(paste(i, " ", sep=""), appendLF=F);
         tg = as.matrix(fread(i));
         t = NULL;
         for(j in 1:dim(genome_size)[1]) {
@@ -81,13 +82,11 @@
         if (length(t)>0) {
             tg = tg[-t,];
         }
-        print(c(dim(g),dim(tg)));
         g = rbind(g, tg);
     }
-    message("Done");
     uchr = sort(unique(as.character(g[,2])));
     g1 = NULL;
-    for(i in uchr){
+    for(i in uchr) {
         t = which(g[,2]==i);
         g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]);
     }
@@ -106,21 +105,19 @@
     }
     cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
     g1 = NULL;
-    message("Generating tracks");
     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));
     for(i in 1:n) {
         tstate = state[,i];
-        t = c(tt,which(tstate[2:L]!=tstate[2:L-1]));
+        t = c(tt, which(tstate[2:L]!=tstate[2:L-1]));
         t = sort(unique(t));
         t0 = c(0, t) + 1;
         t = c(t, L);
         np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
         x = cbind(np[,1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]);
         write.table(as.matrix(x), paste(outpref, i, "bed1", sep="."), quote=F, row.names=F, col.names=F);
-        print(x[1,]);
         system(paste("bedToBigBed ", outpref, ".", i, ".bed1 ", chrom_len_file, " ", outpref, ".", i, ".bb", sep=""));
         system(paste("rm ", paste(outpref, i, "bed1", sep=".")));
     }
@@ -172,3 +169,4 @@
 write.table(track_db, opt$output_track_db, quote=F, row.names=F, col.names=F);
 write.table(c(paste("genome", opt$build), opt$output_build), paste(tracks_output_dir, "genomes_", opt$hub_id, ".txt", sep=""), quote=F, row.names=F, col.names=F);
 write.table(c(paste("hub", opt$hub_id), paste("shortLabel", opt$hub_id), paste("longLabel", opt$hub_name), paste("genomesFile genomes_", opt$hub_id, ".txt", sep=""), "email yzz2 at psu.edu"), opt$output_hub, quote=F, row.names=F, col.names=F);
+