Mercurial > repos > greg > create_genome_tracks

comparison create_genome_tracks.R @ 0:27d5a1d744fd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 10 Nov 2017 13:19:37 -0500
parents
children 64b7d66acd75
comparison
equal deleted inserted replaced
-1:000000000000 0:27d5a1d744fd
1 #!/usr/bin/env Rscript
2
3 suppressPackageStartupMessages(library("optparse"))
4
5 option_list <- list(
6 make_option(c("-b", "--build"), action="store", dest="build", help="Genome build"),
7 make_option(c("-c", "--chrom_len_file"), action="store", dest="chrom_len_file", help="Chromosome length file"),
8 make_option(c("-d", "--header"), action="store", dest="header", default=NULL, help="Track header"),
9 make_option(c("-e", "--state_name"), action="store", dest="state_name", help="State name"),
10 make_option(c("-h", "--hub_id"), action="store", dest="hub_id", help="Not sure what this is"),
11 make_option(c("-i", "--cell_info"), action="store", dest="cell_info", default=NULL, help="Not sure what this is"),
12 make_option(c("-n", "--hub_name"), action="store", dest="hub_name", default=NULL, help="Not sure what this is"),
13 make_option(c("-p", "--parameter_files"), action="store", dest="parameter_files", help="Comma-separated list of IDEAS parameter files"),
14 make_option(c("-s", "--state_files"), action="store", dest="state_files", help="Comma-separated list of IDEAS state files"),
15 make_option(c("-t", "--target_url"), action="store", dest="target_url", help="target url for tracks, not sure what it is used for"),
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"),
18 make_option(c("-w", "--output_hub"), action="store", dest="output_hub", help="Output hub file")
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 #ideas_state_files <- strsplit(opt$state_files, ",")
26
27 get_state_color<-function(statemean, markcolor=NULL)
28 {
29 if(length(markcolor)==0) {
30 markcolor=rep("",dim(statemean)[2]);
31 markcolor[order(apply(statemean,2,sd),decreasing=T)]=hsv((1:dim(statemean)[2]-1)/dim(statemean)[2],1,1);
32 markcolor=t(col2rgb(markcolor));
33 }
34 rg=apply(statemean,1,range);
35 mm=NULL;
36 for(i in 1:dim(statemean)[1]) {
37 mm=rbind(mm,(statemean[i,]-rg[1,i])/(rg[2,i]-rg[1,i]+1e-10));
38 }
39 mm = mm^6;
40 mm = mm / (apply(mm, 1, sum) + 1e-10);
41 mycol = mm%*%markcolor;
42 s = apply(statemean, 1, max);
43 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
44 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
45 h[,2] = h[,2] * s;
46 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
47 rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}), h);
48 return(rt);
49 }
50
51 create_track<-function(state_files, chrom_len_file, outpref, state_color, header, state_name) {
52 message("Reading state file: ", appendLF=FALSE);
53 library("data.table");
54 genomesz = read.table(chrom_len_file);
55 g = NULL;
56 for(i in state_files) {
57 message(paste(i, " ", sep=""), appendLF=F);
58 tg = as.matrix(fread(i));
59 t = NULL;
60 for(j in 1:dim(genomesz)[1]) {
61 t = c(t, which(tg[,2]==as.character(genomesz[j, 1]) & as.numeric(tg[,4]) > as.numeric(genomesz[j, 2])));
62 }
63 if(length(t) > 0) {
64 tg = tg[-t,];
65 }
66 t = which(is.na(match(tg[,2], genomesz[,1]))==T);
67 if(length(t)>0) {
68 tg = tg[-t,];
69 }
70 print(c(dim(g),dim(tg)));
71 g = rbind(g, tg);
72 }
73 message("Done");
74 uchr = sort(unique(as.character(g[,2])));
75 g1 = NULL;
76 for(i in uchr){
77 t = which(g[,2]==i);
78 g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]);
79 }
80 g = NULL;
81 chr = as.character(g1[,2]);
82 posst = as.numeric(g1[,3]);
83 posed = as.numeric(g1[,4]);
84 state = as.matrix(g1[,5:(dim(g1)[2]-1)]);
85 if(length(state_name)==0) {
86 state_name=0:max(state);
87 }
88 L = dim(g1)[1];
89 n = dim(state)[2];
90 if(length(header) > 0) {
91 colnames(g1) = header;
92 }
93 cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]);
94 g1 = NULL;
95 message("Generating tracks");
96 options(scipen=999);
97 tt = which(chr[2:L]!=chr[2:L-1]);
98 tt = c(tt,which(posst[2:L]!=posed[2:L-1]));
99 tt = sort(unique(tt));
100 for(i in 1:n) {
101 tstate = state[,i];
102 t = c(tt,which(tstate[2:L]!=tstate[2:L-1]));
103 t = sort(unique(t));
104 t0 = c(0, t) + 1;
105 t = c(t, L);
106 np = cbind(chr[t], posst[t0], posed[t], tstate[t]);
107 x = cbind(np[,1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]);
108 write.table(as.matrix(x), paste(outpref, i, "bed1", sep="."), quote=F, row.names=F, col.names=F);
109 print(x[1,]);
110 system(paste("bedToBigBed ", outpref, ".", i, ".bed1 ", chrom_len_file, " ", outpref, ".", i, ".bb", sep=""));
111 system(paste("rm ", paste(outpref, i, "bed1", sep=".")));
112 }
113 return(cells);
114 }
115
116 # cell_info: shortid match with those in states,
117 # long id, cell description, text color
118 #run<-function(state_files, hub_id, opt$build, chrom_len_file, state_color, opt$target_url="http://bx.psu.edu/~yuzhang/tmp/", tracks_output_dir=NULL, hub_name=NULL, cell_info=NULL, header=NULL, statename=NULL)
119
120 if(length(opt$hub_name) == 0) {
121 hub_name = opt$hub_id;
122 }
123
124 # Create the tracks output directory.
125 if(length(tracks_output_dir) == 0) {
126 tracks_output_dir = paste("tracks_", hub_id, "/", sep="");
127 }
128 dir.create(tracks_output_dir, showWarnings=FALSE);
129
130 # Create the color scheme.
131 mc = NULL;
132 x = read.table(opt$parameter_files, comment="!", nrows=1);
133 l = as.integer(regexpr("\\*", as.matrix(x)));
134 l = min(which(l>0))-2;
135 x = as.matrix(read.table(opt$parameter_files));
136 if(length(opt$parameter_files) > 1) {
137 for(i in 2:length(opt$parameter_files)) {
138 x=x+as.matrix(read.table(paste(opt$parameter_files[i],".para",sep="")));
139 }
140 }
141 p = x[,1] / sum(x[,1]);
142 m = array(as.matrix(x[,1:l+1] / x[,1]), dim=c(dim(x)[1], l));
143 state_color = get_state_color(m, mc);
144
145 # Create the tracks.
146 cells = create_track(opt$state_files, opt$chrom_len_file, paste(tracks_output_dir, hub_id, sep=""), state_color, header=opt$header, state_name=opt$state_name);
147 cell_info = opt$cell_info
148 if(length(cell_info) == 0) {
149 cell_info = cbind(cells, cells, cells, "#000000");
150 cell_info = array(cell_info, dim=c(length(cells), 4));
151 }
152 cell_info = as.matrix(cell_info);
153 track_db = NULL;
154
155 for(i in 1:length(cells)) {
156 ii = which(cells[i] == cell_info[,1]);
157 if(length(ii) == 0) {
158 next;
159 }
160 ii = ii[1];
161 track_db = c(track_db, paste("track bigBed", i, sep=""));
162 track_db = c(track_db, paste("priority", ii));
163 track_db = c(track_db, "type bigBed 9 .");
164 track_db = c(track_db, "itemRgb on");
165 track_db = c(track_db, "maxItems 100000");
166 track_db = c(track_db, paste("bigDataUrl ", opt$target_url, hub_id, ".", i, ".bb", sep=""));
167 track_db = c(track_db, paste("shortLabel", cell_info[ii, 2]));
168 track_db = c(track_db, paste("longLabel", paste(opt$hub_name, cell_info[ii, 3])));
169 track_db = c(track_db, paste("color", paste(c(col2rgb(cell_info[ii, 4])), collapse=",")));
170 track_db = c(track_db, "visibility dense");
171 track_db = c(track_db, "");
172 }
173
174 # Write the outputs.
175 write.table(track_db, opt$output_track_db, quote=F, row.names=F, col.names=F);
176 write.table(c(paste("genome", opt$build), opt$output_build, paste(tracks_output_dir, "genomes_", hub_id, ".txt", sep=""), quote=F,row.names=F,col.names=F);
177 write.table(c(paste("hub", hub_id), paste("shortLabel", hub_id), paste("longLabel", opt$hub_name), paste("genomesFile genomes_", hub_id, ".txt", sep=""), "email yzz2 at psu.edu"), opt$output_hub, quote=F, row.names=F, col.names=F);