|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
5
|
3 suppressPackageStartupMessages(library("data.table"))
|
|
0
|
4 suppressPackageStartupMessages(library("optparse"))
|
|
|
5
|
|
|
6 option_list <- list(
|
|
|
7 make_option(c("-b", "--build"), action="store", dest="build", help="Genome build"),
|
|
|
8 make_option(c("-c", "--chrom_len_file"), action="store", dest="chrom_len_file", help="Chromosome length file"),
|
|
20
|
9 make_option(c("-e", "--email"), action="store", dest="email", help="User email address"),
|
|
|
10 make_option(c("-l", "--long_label"), action="store", dest="long_label", help="Hub long label"),
|
|
|
11 make_option(c("-n", "--hub_name"), action="store", dest="hub_name", default=NULL, help="Hub name without spaces"),
|
|
0
|
12 make_option(c("-p", "--input_dir_para"), action="store", dest="input_dir_para", help="Directory containing .para outputs from IDEAS"),
|
|
|
13 make_option(c("-q", "--input_dir_state"), action="store", dest="input_dir_state", help="Directory containing .state outputs from IDEAS"),
|
|
20
|
14 make_option(c("-s", "--short_label"), action="store", dest="short_label", help="Hub short label"),
|
|
25
|
15 make_option(c("-t", "--track_color"), action="store", dest="track_color", default=NULL, help="User specified track color"),
|
|
0
|
16 make_option(c("-u", "--output_track_db"), action="store", dest="output_track_db", help="Output track db file"),
|
|
11
|
17 make_option(c("-w", "--output_trackhub"), action="store", dest="output_trackhub", help="Output hub file"),
|
|
|
18 make_option(c("-x", "--output_trackhub_files_path"), action="store", dest="output_trackhub_files_path", help="Output hub extra files path")
|
|
0
|
19 )
|
|
|
20
|
|
|
21 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
22 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
23 opt <- args$options
|
|
|
24
|
|
|
25 create_color_scheme = function(input_dir_para) {
|
|
|
26 # Create the color scheme.
|
|
5
|
27 para_files <- list.files(path=input_dir_para, full.names=TRUE);
|
|
0
|
28 mc = NULL;
|
|
5
|
29 x = read.table(para_files[1], comment="!", nrows=1);
|
|
0
|
30 l = as.integer(regexpr("\\*", as.matrix(x)));
|
|
|
31 l = min(which(l>0)) - 2;
|
|
5
|
32 x = as.matrix(read.table(para_files[1]));
|
|
0
|
33 if (length(para_files) > 1) {
|
|
|
34 for (i in 2:length(para_files)) {
|
|
5
|
35 x = x + as.matrix(read.table(para_files[i]));
|
|
0
|
36 }
|
|
|
37 }
|
|
|
38 p = x[,1] / sum(x[,1]);
|
|
|
39 m = array(as.matrix(x[,1:l+1] / x[,1]), dim=c(dim(x)[1], l));
|
|
25
|
40 track_color <- get_track_color(m, mc);
|
|
|
41 return(track_color);
|
|
0
|
42 }
|
|
|
43
|
|
22
|
44 get_hue = function(val, min=0, max=1) {
|
|
|
45 if (is.null(val)) {
|
|
|
46 val = max;
|
|
|
47 } else if (is.na(val)) {
|
|
|
48 val = max;
|
|
|
49 } else if (val < min) {
|
|
|
50 val = min;
|
|
|
51 } else if (val > max) {
|
|
|
52 val = max;
|
|
|
53 }
|
|
|
54 return(val);
|
|
|
55 }
|
|
|
56
|
|
20
|
57 get_track_file_name = function(base_track_file_name, index, ext) {
|
|
|
58 track_file_name <- paste(base_track_file_name, index, ext, sep=".");
|
|
|
59 return(track_file_name);
|
|
|
60 }
|
|
|
61
|
|
33
|
62 get_track_color = function(statemean, markcolor=NULL) {
|
|
22
|
63 sm_width = dim(statemean)[1];
|
|
|
64 sm_height = dim(statemean)[2];
|
|
0
|
65 if (length(markcolor) == 0) {
|
|
22
|
66 markcolor = rep("", sm_height);
|
|
|
67 markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:sm_height-1)/sm_height, 1, 1);
|
|
0
|
68 markcolor = t(col2rgb(markcolor));
|
|
|
69 }
|
|
|
70 rg = apply(statemean, 1, range);
|
|
|
71 mm = NULL;
|
|
22
|
72 for (i in 1:sm_width) {
|
|
|
73 mm = rbind(mm, statemean[i,] - rg[1,i] + 1e-10 / rg[2,i] - rg[1,i] + 1e-10);
|
|
0
|
74 }
|
|
20
|
75 mm = mm^6;
|
|
|
76 if(dim(mm)[2] > 1) {
|
|
|
77 mm = mm / (apply(mm, 1, sum) + 1e-10);
|
|
|
78 }
|
|
0
|
79 mycol = mm%*%markcolor;
|
|
|
80 s = apply(statemean, 1, max);
|
|
|
81 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
|
|
22
|
82 h = t(apply(mycol, 1, function(x) {rgb2hsv(r=get_hue(x[1]), g=get_hue(x[2]), b=get_hue(x[3]))}));
|
|
0
|
83 h[,2] = h[,2] * s;
|
|
5
|
84 h = apply(h, 1, function(x) {hsv(x[1], x[2], x[3])});
|
|
|
85 rt = cbind(apply(t(col2rgb(h)), 1, function(x) {paste(x, collapse=",")}), h);
|
|
0
|
86 return(rt);
|
|
|
87 }
|
|
|
88
|
|
20
|
89 create_primary_html = function(output_trackhub, tracks_dir) {
|
|
22
|
90 track_files <- list.files(path=tracks_dir);
|
|
29
|
91 s <- paste('<html>\n <head>\n </head>\n <body>\n <ul>\n', sep="");
|
|
20
|
92 for (i in 1:length(track_files)) {
|
|
33
|
93 s <- paste(s, ' <li><a href="', track_files[i], '">', track_files[i], '</a></li>\n', sep="");
|
|
22
|
94 }
|
|
29
|
95 s <- paste(s, ' </ul>\n </body>\n</html>', sep="");
|
|
22
|
96 cat(s, file=output_trackhub);
|
|
20
|
97 }
|
|
|
98
|
|
25
|
99 create_track = function(input_dir_state, chrom_len_file, base_track_file_name, track_color) {
|
|
5
|
100 state_files <- list.files(path=input_dir_state, full.names=TRUE);
|
|
0
|
101 genome_size = read.table(chrom_len_file);
|
|
|
102 g = NULL;
|
|
|
103 for(i in state_files) {
|
|
|
104 tg = as.matrix(fread(i));
|
|
|
105 t = NULL;
|
|
|
106 for(j in 1:dim(genome_size)[1]) {
|
|
|
107 t = c(t, which(tg[,2]==as.character(genome_size[j, 1]) & as.numeric(tg[,4]) > as.numeric(genome_size[j, 2])));
|
|
33
|
108 }
|
|
0
|
109 if (length(t) > 0) {
|
|
|
110 tg = tg[-t,];
|
|
|
111 }
|
|
|
112 t = which(is.na(match(tg[,2], genome_size[,1]))==T);
|
|
|
113 if (length(t)>0) {
|
|
|
114 tg = tg[-t,];
|
|
|
115 }
|
|
|
116 g = rbind(g, tg);
|
|
|
117 }
|
|
|
118 uchr = sort(unique(as.character(g[,2])));
|
|
|
119 g1 = NULL;
|
|
5
|
120 for(i in uchr) {
|
|
0
|
121 t = which(g[,2]==i);
|
|
|
122 g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]);
|
|
|
123 }
|
|
|
124 g = NULL;
|
|
|
125 chr = as.character(g1[,2]);
|
|
|
126 posst = as.numeric(g1[,3]);
|
|
|
127 posed = as.numeric(g1[,4]);
|
|
|
128 state = as.matrix(g1[,5:(dim(g1)[2]-1)]);
|
|
20
|
129 state_name = 0:max(state);
|
|
0
|
130 L = dim(g1)[1];
|
|
|
131 n = dim(state)[2];
|
|
|
132 cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
|
|
|
133 g1 = NULL;
|
|
|
134 options(scipen=999);
|
|
|
135 tt = which(chr[2:L]!=chr[2:L-1]);
|
|
11
|
136 tt = c(tt, which(posst[2:L]!=posed[2:L-1]));
|
|
0
|
137 tt = sort(unique(tt));
|
|
|
138 for(i in 1:n) {
|
|
|
139 tstate = state[,i];
|
|
5
|
140 t = c(tt, which(tstate[2:L]!=tstate[2:L-1]));
|
|
0
|
141 t = sort(unique(t));
|
|
|
142 t0 = c(0, t) + 1;
|
|
|
143 t = c(t, L);
|
|
|
144 np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
|
|
25
|
145 x = cbind(np[,1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], track_color[as.numeric(np[,4])+1]);
|
|
20
|
146 track_file_name_bed <- get_track_file_name(base_track_file_name, i, "bed");
|
|
|
147 track_file_name_bigbed <- get_track_file_name(base_track_file_name, i, "bigbed");
|
|
11
|
148 write.table(as.matrix(x), track_file_name_bed, quote=F, row.names=F, col.names=F);
|
|
|
149 system(paste("bedToBigBed ", track_file_name_bed, chrom_len_file, " ", track_file_name_bigbed));
|
|
15
|
150 system(paste("rm ", track_file_name_bed));
|
|
0
|
151 }
|
|
|
152 return(cells);
|
|
|
153 }
|
|
|
154
|
|
25
|
155 create_track_db = function(input_dir_state, chrom_len_file, tracks_dir, hub_name, short_label, long_label, track_color) {
|
|
20
|
156 base_track_file_name <- paste(tracks_dir, hub_name, sep="");
|
|
25
|
157 cells = create_track(input_dir_state, chrom_len_file, base_track_file_name, track_color);
|
|
20
|
158 cell_info = cbind(cells, cells, cells, "#000000");
|
|
|
159 cell_info = array(cell_info, dim=c(length(cells), 4));
|
|
0
|
160 cell_info = as.matrix(cell_info);
|
|
|
161 track_db = NULL;
|
|
25
|
162 for (i in 1:length(cells)) {
|
|
20
|
163
|
|
22
|
164 ii = which(cells[i] == cell_info[,1]);
|
|
0
|
165 if (length(ii) == 0) {
|
|
|
166 next;
|
|
|
167 }
|
|
|
168 ii = ii[1];
|
|
20
|
169 # trackDb.txt track entry.
|
|
32
|
170 track_db = c(track_db, paste("hub ", hub_name, "_track_", i, sep=""));
|
|
20
|
171 track_db = c(track_db, "type bigBed");
|
|
22
|
172 track_db = c(track_db, paste("bigDataUrl", get_track_file_name(base_track_file_name, i, "bigbed"), sep=" "));
|
|
|
173 track_db = c(track_db, paste("shortLabel", short_label, sep=" "));
|
|
|
174 track_db = c(track_db, paste("longLabel", long_label, sep=" "));
|
|
|
175 track_db = c(track_db, paste("priority", ii));
|
|
0
|
176 track_db = c(track_db, "itemRgb on");
|
|
|
177 track_db = c(track_db, "maxItems 100000");
|
|
22
|
178 track_db = c(track_db, paste("color", paste(c(col2rgb(cell_info[ii,4])), collapse=","), sep=" "));
|
|
0
|
179 track_db = c(track_db, "visibility dense");
|
|
33
|
180 track_db = c(track_db, "");
|
|
0
|
181 }
|
|
11
|
182 return(track_db);
|
|
0
|
183 }
|
|
|
184
|
|
|
185 # Create the color scheme.
|
|
28
|
186 if (length(opt$track_color) == 0) {
|
|
25
|
187 track_color <- create_color_scheme(opt$input_dir_para);
|
|
28
|
188 } else {
|
|
|
189 track_color <- col2rgb(opt$track_color);
|
|
25
|
190 }
|
|
0
|
191
|
|
20
|
192 # Create the hub.txt output.
|
|
34
|
193 hub_name_line <- c(paste("hub ", opt$hub_name), sep = "");
|
|
|
194 short_label_line <- c(paste("shortLabel ", opt$short_label), sep = "");
|
|
|
195 long_label_line <- c(paste("longLabel ", opt$longLabel), sep = "");
|
|
|
196 genomes_txt_line <- c(paste("genomesFile genomes.txt"), sep="");
|
|
|
197 email_line <- c(paste("email ", opt$email), sep = "");
|
|
33
|
198 contents <- c(paste(hub_name_line, short_label_line, long_label_line, genomes_txt_line, email_line), sep="");
|
|
34
|
199 hub_dir <- c(paste(opt$output_trackhub_files_path, "/", "myHub", "/"), sep="");
|
|
22
|
200 dir.create(hub_dir, showWarnings=FALSE);
|
|
34
|
201 hub_file_path <- c(paste(hub_dir, "hub.txt"), sep="");
|
|
20
|
202 write.table(contents, file=hub_file_path, quote=F, row.names=F, col.names=F);
|
|
|
203
|
|
|
204 # Create the genomes.txt output.
|
|
34
|
205 genome_line <- c(paste("genome", opt$build), sep="");
|
|
|
206 track_db_line <- c(paste("trackDb ", opt$build, "/", "trackDb.txt"), sep="");
|
|
33
|
207 contents <- c(paste(genome_line, track_db_line), sep="");
|
|
30
|
208 genomes_file_path <- paste(hub_dir, "genomes.txt", sep="");
|
|
20
|
209 write.table(contents, file=genomes_file_path, quote=F, row.names=F, col.names=F);
|
|
|
210
|
|
0
|
211 # Create the tracks.
|
|
34
|
212 tracks_dir <- c(paste(hub_dir, opt$build, "/"), sep="");
|
|
11
|
213 dir.create(tracks_dir, showWarnings=FALSE);
|
|
25
|
214 track_db <- create_track_db(opt$input_dir_state, opt$chrom_len_file, tracks_dir, opt$hub_name, opt$short_label, opt$long_label, track_color);
|
|
0
|
215
|
|
20
|
216 # Create the trackDb.txt output.
|
|
34
|
217 track_db_file_path <- c(paste(tracks_dir, "trackDb.txt"), sep="");
|
|
13
|
218 write.table(track_db, file=track_db_file_path, quote=F, row.names=F, col.names=F);
|
|
11
|
219
|
|
20
|
220 # Create the primary HTML dataset.
|
|
|
221 create_primary_html(opt$output_trackhub, tracks_dir);
|
