|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
5
|
3 suppressPackageStartupMessages(library("data.table"))
|
|
0
|
4 suppressPackageStartupMessages(library("optparse"))
|
|
44
|
5 suppressPackageStartupMessages(library("viridislite"))
|
|
0
|
6
|
|
|
7 option_list <- list(
|
|
39
|
8 make_option(c("--build"), action="store", dest="build", help="Genome build"),
|
|
|
9 make_option(c("--chrom_len_file"), action="store", dest="chrom_len_file", help="Chromosome length file"),
|
|
|
10 make_option(c("--email"), action="store", dest="email", help="User email address"),
|
|
|
11 make_option(c("--galaxy_url"), action="store", dest="galaxy_url", help="Galaxy instance base URL"),
|
|
|
12 make_option(c("--hub_name"), action="store", dest="hub_name", default=NULL, help="Hub name without spaces"),
|
|
|
13 make_option(c("--input_dir_para"), action="store", dest="input_dir_para", help="Directory containing .para outputs from IDEAS"),
|
|
|
14 make_option(c("--input_dir_state"), action="store", dest="input_dir_state", help="Directory containing .state outputs from IDEAS"),
|
|
|
15 make_option(c("--long_label"), action="store", dest="long_label", help="Hub long label"),
|
|
|
16 make_option(c("--output_trackhub"), action="store", dest="output_trackhub", help="Output hub file"),
|
|
|
17 make_option(c("--output_trackhub_files_path"), action="store", dest="output_trackhub_files_path", help="Output hub extra files path"),
|
|
|
18 make_option(c("--output_trackhub_id"), action="store", dest="output_trackhub_id", help="Encoded output_trackhub dataset id"),
|
|
42
|
19 make_option(c("--short_label"), action="store", dest="short_label", help="Hub short label")
|
|
0
|
20 )
|
|
|
21
|
|
|
22 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
23 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
24 opt <- args$options
|
|
|
25
|
|
39
|
26 create_primary_html = function(output_trackhub, tracks_dir, build) {
|
|
|
27 track_files <- list.files(path=tracks_dir);
|
|
|
28 s <- paste('<html><head></head><body>', sep="\n");
|
|
|
29 s <- paste(s, '<h3>Contents of directory ~/myHub/', build, ' required by UCSC TrackHub</h3>\n', sep="");
|
|
|
30 s <- paste(s, '<ul>\n', sep="")
|
|
|
31 for (i in 1:length(track_files)) {
|
|
|
32 s <- paste(s, '<li><a href="', 'myHub/', build, "/", track_files[i], '">', track_files[i], '</a></li>\n', sep="");
|
|
0
|
33 }
|
|
39
|
34 s <- paste(s, '</ul>\n</body>\n</html>', sep="");
|
|
22
|
35 cat(s, file=output_trackhub);
|
|
20
|
36 }
|
|
|
37
|
|
42
|
38 create_track = function(input_dir_para, input_dir_state, chrom_len_file, base_track_file_name) {
|
|
|
39 para_files <- list.files(path=input_dir_para, full.names=TRUE);
|
|
5
|
40 state_files <- list.files(path=input_dir_state, full.names=TRUE);
|
|
0
|
41 genome_size = read.table(chrom_len_file);
|
|
|
42 g = NULL;
|
|
|
43 for(i in state_files) {
|
|
|
44 tg = as.matrix(fread(i));
|
|
|
45 t = NULL;
|
|
|
46 for(j in 1:dim(genome_size)[1]) {
|
|
|
47 t = c(t, which(tg[,2]==as.character(genome_size[j, 1]) & as.numeric(tg[,4]) > as.numeric(genome_size[j, 2])));
|
|
33
|
48 }
|
|
0
|
49 if (length(t) > 0) {
|
|
|
50 tg = tg[-t,];
|
|
|
51 }
|
|
|
52 t = which(is.na(match(tg[,2], genome_size[,1]))==T);
|
|
|
53 if (length(t)>0) {
|
|
|
54 tg = tg[-t,];
|
|
|
55 }
|
|
|
56 g = rbind(g, tg);
|
|
|
57 }
|
|
|
58 uchr = sort(unique(as.character(g[,2])));
|
|
|
59 g1 = NULL;
|
|
5
|
60 for(i in uchr) {
|
|
0
|
61 t = which(g[,2]==i);
|
|
|
62 g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]);
|
|
|
63 }
|
|
|
64 g = NULL;
|
|
|
65 chr = as.character(g1[,2]);
|
|
|
66 posst = as.numeric(g1[,3]);
|
|
|
67 posed = as.numeric(g1[,4]);
|
|
|
68 state = as.matrix(g1[,5:(dim(g1)[2]-1)]);
|
|
20
|
69 state_name = 0:max(state);
|
|
0
|
70 L = dim(g1)[1];
|
|
|
71 n = dim(state)[2];
|
|
|
72 cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
|
|
|
73 g1 = NULL;
|
|
|
74 options(scipen=999);
|
|
|
75 tt = which(chr[2:L]!=chr[2:L-1]);
|
|
11
|
76 tt = c(tt, which(posst[2:L]!=posed[2:L-1]));
|
|
0
|
77 tt = sort(unique(tt));
|
|
|
78 for(i in 1:n) {
|
|
42
|
79 state_color <- get_state_color(para_files[i])
|
|
0
|
80 tstate = state[,i];
|
|
5
|
81 t = c(tt, which(tstate[2:L]!=tstate[2:L-1]));
|
|
0
|
82 t = sort(unique(t));
|
|
|
83 t0 = c(0, t) + 1;
|
|
|
84 t = c(t, L);
|
|
|
85 np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
|
|
20
|
86 track_file_name_bed <- get_track_file_name(base_track_file_name, i, "bed");
|
|
|
87 track_file_name_bigbed <- get_track_file_name(base_track_file_name, i, "bigbed");
|
|
42
|
88 x = cbind(np[, 1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]);
|
|
11
|
89 write.table(as.matrix(x), track_file_name_bed, quote=F, row.names=F, col.names=F);
|
|
|
90 system(paste("bedToBigBed ", track_file_name_bed, chrom_len_file, " ", track_file_name_bigbed));
|
|
15
|
91 system(paste("rm ", track_file_name_bed));
|
|
0
|
92 }
|
|
|
93 return(cells);
|
|
|
94 }
|
|
|
95
|
|
42
|
96 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) {
|
|
20
|
97 base_track_file_name <- paste(tracks_dir, hub_name, sep="");
|
|
42
|
98 cells = create_track(input_dir_para, input_dir_state, chrom_len_file, base_track_file_name);
|
|
20
|
99 cell_info = cbind(cells, cells, cells, "#000000");
|
|
|
100 cell_info = array(cell_info, dim=c(length(cells), 4));
|
|
0
|
101 cell_info = as.matrix(cell_info);
|
|
|
102 track_db = NULL;
|
|
25
|
103 for (i in 1:length(cells)) {
|
|
39
|
104 big_data_url <- get_big_data_url(galaxy_url, encoded_dataset_id, tracks_dir, i, build);
|
|
22
|
105 ii = which(cells[i] == cell_info[,1]);
|
|
0
|
106 if (length(ii) == 0) {
|
|
|
107 next;
|
|
|
108 }
|
|
|
109 ii = ii[1];
|
|
20
|
110 # trackDb.txt track entry.
|
|
32
|
111 track_db = c(track_db, paste("hub ", hub_name, "_track_", i, sep=""));
|
|
20
|
112 track_db = c(track_db, "type bigBed");
|
|
39
|
113 track_db = c(track_db, paste("bigDataUrl", big_data_url, sep=" "));
|
|
22
|
114 track_db = c(track_db, paste("shortLabel", short_label, sep=" "));
|
|
|
115 track_db = c(track_db, paste("longLabel", long_label, sep=" "));
|
|
|
116 track_db = c(track_db, paste("priority", ii));
|
|
0
|
117 track_db = c(track_db, "itemRgb on");
|
|
|
118 track_db = c(track_db, "maxItems 100000");
|
|
22
|
119 track_db = c(track_db, paste("color", paste(c(col2rgb(cell_info[ii,4])), collapse=","), sep=" "));
|
|
0
|
120 track_db = c(track_db, "visibility dense");
|
|
33
|
121 track_db = c(track_db, "");
|
|
0
|
122 }
|
|
11
|
123 return(track_db);
|
|
0
|
124 }
|
|
|
125
|
|
39
|
126 get_big_data_url = function(galaxy_url, encoded_dataset_id, tracks_dir, index, build) {
|
|
|
127 track_files <- list.files(path=tracks_dir, pattern="\\.bigbed");
|
|
|
128 s <- paste(galaxy_url, 'datasets/', encoded_dataset_id, '/display/myHub/', build, '/', track_files[index], sep="");
|
|
|
129 return(s)
|
|
|
130 }
|
|
|
131
|
|
42
|
132 get_rgb<-function(statemean, markcolor=NULL)
|
|
44
|
133 {
|
|
42
|
134 if(length(markcolor) == 0) {
|
|
|
135 markcolor = rep("",dim(statemean)[2]);
|
|
|
136 markcolor[order(apply(statemean,2,sd),decreasing=T)]=hsv((1:dim(statemean)[2]-1)/dim(statemean)[2],1,1)
|
|
39
|
137 markcolor = t(col2rgb(markcolor));
|
|
|
138 }
|
|
44
|
139
|
|
39
|
140 rg = apply(statemean, 1, range);
|
|
|
141 mm = NULL;
|
|
42
|
142 for(i in 1:dim(statemean)[1]) {
|
|
|
143 mm = rbind(mm, (statemean[i,]-rg[1, i]+1e-10)/(rg[2, i]-rg[1, i]+1e-10));
|
|
39
|
144 }
|
|
44
|
145 mm = mm^5;
|
|
42
|
146 if (dim(mm)[2]>1) {
|
|
|
147 mm = mm / (apply(mm, 1, sum)+1e-10);
|
|
39
|
148 }
|
|
|
149 mycol = mm%*%markcolor;
|
|
|
150 s = apply(statemean, 1, max);
|
|
42
|
151 s = (s-min(s))/(max(s)-min(s)+1e-10);
|
|
|
152 mycol = round(255-(255-mycol)*s/0.5);
|
|
|
153 mycol[mycol<0] = 0;
|
|
|
154 rt = paste(mycol[,1], mycol[,2], mycol[,3], sep=",");
|
|
|
155 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
|
|
|
156 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
|
|
|
157 rt = cbind(rt, h);
|
|
|
158 return(rt);
|
|
44
|
159
|
|
42
|
160 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
|
|
|
161 h[,2] = h[,2]*s;
|
|
|
162 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
|
|
|
163 rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}), h);
|
|
39
|
164 return(rt);
|
|
|
165 }
|
|
|
166
|
|
42
|
167 get_state_color <- function(para_file, cols=c("white", "dark blue")) {
|
|
|
168 x = read.table(parafile, comment.char="#", header=T);
|
|
|
169 k = dim(x)[2];
|
|
|
170 l = dim(x)[1];
|
|
|
171 p = (sqrt(9+8*(k-1))-3)/2;
|
|
|
172 m = as.matrix(x[,1+1:p]/x[,1]);
|
|
|
173 colnames(m) = colnames(x)[1+1:p];
|
|
|
174 marks = colnames(m);
|
|
|
175 rownames(m) = paste(1:l-1," (",round(x[,1]/sum(x[,1])*10000)/100,"%)",sep="");
|
|
|
176 par(mar=c(6,1,1,6));
|
|
|
177 rg = range(m);
|
|
|
178 colors = 0:100/100*(rg[2]-rg[1])+rg[1];
|
|
|
179 my_palette = colorRampPalette(cols)(n=100);
|
|
|
180 defpalette = palette(my_palette);
|
|
|
181 plot(NA,NA,xlim=c(0,p+0.7),ylim=c(0,l),xaxt="n",yaxt="n",xlab=NA,ylab=NA,frame.plot=F);
|
|
|
182 axis(1,at=1:p-0.5,labels=colnames(m),las=2);
|
|
|
183 axis(4,at=1:l-0.5,labels=rownames(m),las=2);
|
|
|
184 rect(rep(1:p-1,l),rep(1:l-1,each=p),rep(1:p,l),rep(1:l,each=p),col=round((t(m)-rg[1])/(rg[2]-rg[1])*100));
|
|
44
|
185 markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
|
|
42
|
186 for(i in 1:length(marks)) {
|
|
|
187 if (regexpr("h3k4me3",tolower(marks[i]))>0) {
|
|
|
188 markcolor[i,]=c(255,0,0);
|
|
|
189 }
|
|
|
190 if (regexpr("h3k4me2",tolower(marks[i]))>0) {
|
|
|
191 markcolor[i,]=c(250,100,0);
|
|
|
192 }
|
|
|
193 if (regexpr("h3k4me1",tolower(marks[i]))>0) {
|
|
|
194 markcolor[i,]=c(250,250,0);
|
|
|
195 }
|
|
|
196 if (regexpr("h3k36me3",tolower(marks[i]))>0) {
|
|
|
197 markcolor[i,]=c(0,150,0);
|
|
|
198 }
|
|
|
199 if (regexpr("h2a",tolower(marks[i]))>0) {
|
|
|
200 markcolor[i,]=c(0,150,150);
|
|
|
201 }
|
|
|
202 if (regexpr("dnase",tolower(marks[i]))>0) {
|
|
|
203 markcolor[i,]=c(0,200,200);
|
|
|
204 }
|
|
|
205 if (regexpr("atac",tolower(marks[i]))>0) {
|
|
|
206 markcolor[i,]=c(0,200,200);
|
|
|
207 }
|
|
|
208 if (regexpr("h3k9ac",tolower(marks[i]))>0) {
|
|
|
209 markcolor[i,]=c(250,0,200);
|
|
|
210 }
|
|
|
211 if (regexpr("h3k9me3",tolower(marks[i]))>0) {
|
|
|
212 markcolor[i,]=c(100,100,100);
|
|
|
213 }
|
|
|
214 if (regexpr("h3k27ac",tolower(marks[i]))>0) {
|
|
|
215 markcolor[i,]=c(250,150,0);
|
|
|
216 }
|
|
|
217 if (regexpr("h3k27me3",tolower(marks[i]))>0) {
|
|
|
218 markcolor[i,]=c(0,0,225);
|
|
|
219 }
|
|
|
220 if (regexpr("h3k79me2",tolower(marks[i]))>0) {
|
|
|
221 markcolor[i,]=c(200,0,200);
|
|
|
222 }
|
|
|
223 if (regexpr("h4k20me1",tolower(marks[i]))>0) {
|
|
|
224 markcolor[i,]=c(50,200,50);
|
|
|
225 }
|
|
|
226 if (regexpr("ctcf",tolower(marks[i]))>0) {
|
|
|
227 markcolor[i,]=c(200,0,250);
|
|
|
228 }
|
|
|
229 }
|
|
|
230 statecolor = get_rgb(m, markcolor)[,];
|
|
|
231 rect(rep(p+0.2,l),1:l-0.8,rep(p+0.8,l),1:l-0.2,col=statecolor[,2]);
|
|
|
232 palette(defpalette);
|
|
|
233 return(statecolor);
|
|
|
234 }
|
|
|
235
|
|
39
|
236 get_track_file_name = function(base_track_file_name, index, ext) {
|
|
|
237 track_file_name <- paste(base_track_file_name, index, ext, sep=".");
|
|
|
238 return(track_file_name);
|
|
|
239 }
|
|
|
240
|
|
20
|
241 # Create the hub.txt output.
|
|
39
|
242 hub_name_line <- paste("hub ", opt$hub_name, sep="");
|
|
|
243 short_label_line <- paste("shortLabel ", opt$short_label, sep="");
|
|
|
244 long_label_line <- paste("longLabel ", opt$long_label, sep="");
|
|
35
|
245 genomes_txt_line <- paste("genomesFile genomes.txt", sep="");
|
|
39
|
246 email_line <- paste("email ", opt$email, sep="");
|
|
|
247 contents <- paste(hub_name_line, short_label_line, long_label_line, genomes_txt_line, email_line, sep="\n");
|
|
35
|
248 hub_dir <- paste(opt$output_trackhub_files_path, "/", "myHub", "/", sep="");
|
|
22
|
249 dir.create(hub_dir, showWarnings=FALSE);
|
|
35
|
250 hub_file_path <- paste(hub_dir, "hub.txt", sep="");
|
|
20
|
251 write.table(contents, file=hub_file_path, quote=F, row.names=F, col.names=F);
|
|
|
252
|
|
|
253 # Create the genomes.txt output.
|
|
39
|
254 genome_line <- paste("genome ", opt$build, sep="");
|
|
35
|
255 track_db_line <- paste("trackDb ", opt$build, "/", "trackDb.txt", sep="");
|
|
39
|
256 contents <- paste(genome_line, track_db_line, sep="\n");
|
|
30
|
257 genomes_file_path <- paste(hub_dir, "genomes.txt", sep="");
|
|
20
|
258 write.table(contents, file=genomes_file_path, quote=F, row.names=F, col.names=F);
|
|
|
259
|
|
0
|
260 # Create the tracks.
|
|
36
|
261 tracks_dir <- paste(hub_dir, opt$build, "/", sep="");
|
|
11
|
262 dir.create(tracks_dir, showWarnings=FALSE);
|
|
42
|
263 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);
|
|
0
|
264
|
|
20
|
265 # Create the trackDb.txt output.
|
|
35
|
266 track_db_file_path <- paste(tracks_dir, "trackDb.txt", sep="");
|
|
13
|
267 write.table(track_db, file=track_db_file_path, quote=F, row.names=F, col.names=F);
|
|
11
|
268
|
|
20
|
269 # Create the primary HTML dataset.
|
|
42
|
270 create_primary_html(opt$output_trackhub, tracks_dir, opt$build);
|
