|
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"),
|
|
|
9 make_option(c("-d", "--header"), action="store", dest="header", default=NULL, help="Track header"),
|
|
|
10 make_option(c("-e", "--state_name"), action="store", dest="state_name", help="State name"),
|
|
|
11 make_option(c("-f", "--hub_id"), action="store", dest="hub_id", help="Not sure what this is"),
|
|
|
12 make_option(c("-i", "--cell_info"), action="store", dest="cell_info", default=NULL, help="Not sure what this is"),
|
|
|
13 make_option(c("-n", "--hub_name"), action="store", dest="hub_name", default=NULL, help="Not sure what this is"),
|
|
|
14 make_option(c("-p", "--input_dir_para"), action="store", dest="input_dir_para", help="Directory containing .para outputs from IDEAS"),
|
|
|
15 make_option(c("-q", "--input_dir_state"), action="store", dest="input_dir_state", help="Directory containing .state outputs from IDEAS"),
|
|
|
16 make_option(c("-u", "--output_track_db"), action="store", dest="output_track_db", help="Output track db file"),
|
|
|
17 make_option(c("-v", "--output_build"), action="store", dest="output_build", help="Output genome build file"),
|
|
9
|
18 make_option(c("-w", "--output_hub"), action="store", dest="output_hub", help="Output hub file"),
|
|
|
19 make_option(c("-x", "--output_hub_files_path"), action="store", dest="output_hub_files_path", help="Output hub extra files path")
|
|
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
|
|
|
26 create_color_scheme = function(input_dir_para) {
|
|
|
27 # Create the color scheme.
|
|
5
|
28 para_files <- list.files(path=input_dir_para, full.names=TRUE);
|
|
0
|
29 mc = NULL;
|
|
5
|
30 x = read.table(para_files[1], comment="!", nrows=1);
|
|
0
|
31 l = as.integer(regexpr("\\*", as.matrix(x)));
|
|
|
32 l = min(which(l>0)) - 2;
|
|
5
|
33 x = as.matrix(read.table(para_files[1]));
|
|
0
|
34 if (length(para_files) > 1) {
|
|
|
35 for (i in 2:length(para_files)) {
|
|
5
|
36 x = x + as.matrix(read.table(para_files[i]));
|
|
0
|
37 }
|
|
|
38 }
|
|
|
39 p = x[,1] / sum(x[,1]);
|
|
|
40 m = array(as.matrix(x[,1:l+1] / x[,1]), dim=c(dim(x)[1], l));
|
|
5
|
41 state_color <- get_state_color(m, mc);
|
|
|
42 state_color
|
|
0
|
43 }
|
|
|
44
|
|
|
45 get_state_color = function(statemean, markcolor=NULL)
|
|
|
46 {
|
|
|
47 if (length(markcolor) == 0) {
|
|
|
48 markcolor = rep("", dim(statemean)[2]);
|
|
|
49 markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:dim(statemean)[2]-1)/dim(statemean)[2], 1, 1);
|
|
|
50 markcolor = t(col2rgb(markcolor));
|
|
|
51 }
|
|
|
52 rg = apply(statemean, 1, range);
|
|
|
53 mm = NULL;
|
|
|
54 for(i in 1:dim(statemean)[1]) {
|
|
|
55 mm = rbind(mm, (statemean[i,]-rg[1,i])/(rg[2,i]-rg[1,i]+1e-10));
|
|
|
56 }
|
|
|
57 mm = mm^6;
|
|
|
58 mm = mm / (apply(mm, 1, sum) + 1e-10);
|
|
|
59 mycol = mm%*%markcolor;
|
|
|
60 s = apply(statemean, 1, max);
|
|
|
61 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
|
|
5
|
62 h = t(apply(mycol, 1, function(x) {rgb2hsv(x[1], x[2], x[3])}));
|
|
0
|
63 h[,2] = h[,2] * s;
|
|
5
|
64 h = apply(h, 1, function(x) {hsv(x[1], x[2], x[3])});
|
|
|
65 rt = cbind(apply(t(col2rgb(h)), 1, function(x) {paste(x, collapse=",")}), h);
|
|
0
|
66 return(rt);
|
|
|
67 }
|
|
|
68
|
|
|
69 create_track = function(input_dir_state, chrom_len_file, outpref, state_color, header, state_name) {
|
|
5
|
70 state_files <- list.files(path=input_dir_state, full.names=TRUE);
|
|
0
|
71 genome_size = read.table(chrom_len_file);
|
|
|
72 g = NULL;
|
|
|
73 for(i in state_files) {
|
|
|
74 tg = as.matrix(fread(i));
|
|
|
75 t = NULL;
|
|
|
76 for(j in 1:dim(genome_size)[1]) {
|
|
|
77 t = c(t, which(tg[,2]==as.character(genome_size[j, 1]) & as.numeric(tg[,4]) > as.numeric(genome_size[j, 2])));
|
|
|
78 }
|
|
|
79 if (length(t) > 0) {
|
|
|
80 tg = tg[-t,];
|
|
|
81 }
|
|
|
82 t = which(is.na(match(tg[,2], genome_size[,1]))==T);
|
|
|
83 if (length(t)>0) {
|
|
|
84 tg = tg[-t,];
|
|
|
85 }
|
|
|
86 g = rbind(g, tg);
|
|
|
87 }
|
|
|
88 uchr = sort(unique(as.character(g[,2])));
|
|
|
89 g1 = NULL;
|
|
5
|
90 for(i in uchr) {
|
|
0
|
91 t = which(g[,2]==i);
|
|
|
92 g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]);
|
|
|
93 }
|
|
|
94 g = NULL;
|
|
|
95 chr = as.character(g1[,2]);
|
|
|
96 posst = as.numeric(g1[,3]);
|
|
|
97 posed = as.numeric(g1[,4]);
|
|
|
98 state = as.matrix(g1[,5:(dim(g1)[2]-1)]);
|
|
|
99 if (length(state_name)==0) {
|
|
|
100 state_name=0:max(state);
|
|
|
101 }
|
|
|
102 L = dim(g1)[1];
|
|
|
103 n = dim(state)[2];
|
|
|
104 if (length(header) > 0) {
|
|
|
105 colnames(g1) = header;
|
|
|
106 }
|
|
|
107 cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
|
|
|
108 g1 = NULL;
|
|
|
109 options(scipen=999);
|
|
|
110 tt = which(chr[2:L]!=chr[2:L-1]);
|
|
|
111 tt = c(tt,which(posst[2:L]!=posed[2:L-1]));
|
|
|
112 tt = sort(unique(tt));
|
|
|
113 for(i in 1:n) {
|
|
|
114 tstate = state[,i];
|
|
5
|
115 t = c(tt, which(tstate[2:L]!=tstate[2:L-1]));
|
|
0
|
116 t = sort(unique(t));
|
|
|
117 t0 = c(0, t) + 1;
|
|
|
118 t = c(t, L);
|
|
|
119 np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
|
|
|
120 x = cbind(np[,1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]);
|
|
|
121 write.table(as.matrix(x), paste(outpref, i, "bed1", sep="."), quote=F, row.names=F, col.names=F);
|
|
|
122 system(paste("bedToBigBed ", outpref, ".", i, ".bed1 ", chrom_len_file, " ", outpref, ".", i, ".bb", sep=""));
|
|
|
123 system(paste("rm ", paste(outpref, i, "bed1", sep=".")));
|
|
|
124 }
|
|
|
125 return(cells);
|
|
|
126 }
|
|
|
127
|
|
|
128 create_track_db = function(input_dir_state, chrom_len_file, track_file_name, state_color, header, state_name, hub_name, cell_info) {
|
|
|
129 cells = create_track(input_dir_state, chrom_len_file, track_file_name, state_color, header, state_name);
|
|
|
130 if (length(cell_info) == 0) {
|
|
|
131 cell_info = cbind(cells, cells, cells, "#000000");
|
|
|
132 cell_info = array(cell_info, dim=c(length(cells), 4));
|
|
|
133 }
|
|
|
134 cell_info = as.matrix(cell_info);
|
|
|
135 track_db = NULL;
|
|
|
136 for(i in 1:length(cells)) {
|
|
|
137 ii = which(cells[i] == cell_info[,1]);
|
|
|
138 if (length(ii) == 0) {
|
|
|
139 next;
|
|
|
140 }
|
|
|
141 ii = ii[1];
|
|
|
142 track_db = c(track_db, paste("track bigBed", i, sep=""));
|
|
|
143 track_db = c(track_db, paste("priority", ii));
|
|
|
144 track_db = c(track_db, "type bigBed 9 .");
|
|
|
145 track_db = c(track_db, "itemRgb on");
|
|
|
146 track_db = c(track_db, "maxItems 100000");
|
|
|
147 track_db = c(track_db, paste("shortLabel", cell_info[ii, 2]));
|
|
|
148 track_db = c(track_db, paste("longLabel", paste(hub_name, cell_info[ii, 3])));
|
|
|
149 track_db = c(track_db, paste("color", paste(c(col2rgb(cell_info[ii, 4])), collapse=",")));
|
|
|
150 track_db = c(track_db, "visibility dense");
|
|
|
151 track_db = c(track_db, "");
|
|
|
152 }
|
|
|
153 }
|
|
|
154
|
|
|
155 if (length(opt$hub_name) == 0) {
|
|
9
|
156 hub_name <- opt$hub_id;
|
|
|
157 } else {
|
|
|
158 hub_name <- opt$hub_name;
|
|
0
|
159 }
|
|
|
160
|
|
|
161 # Create the tracks output directory.
|
|
|
162 tracks_output_dir = paste("tracks_", opt$hub_id, "/", sep="");
|
|
|
163 dir.create(tracks_output_dir, showWarnings=FALSE);
|
|
|
164
|
|
|
165 # Create the color scheme.
|
|
|
166 state_color <- create_color_scheme(opt$input_dir_para);
|
|
|
167
|
|
|
168 # Create the tracks.
|
|
9
|
169 track_db <- create_track_db(opt$input_dir_state, opt$chrom_len_file, paste(tracks_output_dir, opt$hub_id, sep=""), state_color, opt$header, opt$state_name, hub_name, opt$cell_info);
|
|
0
|
170
|
|
|
171 # Write the outputs.
|
|
|
172 write.table(track_db, opt$output_track_db, quote=F, row.names=F, col.names=F);
|
|
|
173 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);
|
|
9
|
174 write.table(c(paste("hub", opt$hub_id), paste("shortLabel", opt$hub_id), paste("longLabel", hub_name), paste("genomesFile genomes_", opt$hub_id, ".txt", sep=""), "email yzz2 at psu.edu"), paste(opt$output_hub_files_path, "hub_", opt$hub_id, ".txt", sep=""), quote=F, row.names=F, col.names=F);
|
