|
150
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("data.table"))
|
|
|
4 suppressPackageStartupMessages(library("optparse"))
|
|
|
5
|
|
|
6 option_list <- list(
|
|
|
7 make_option(c("--burnin_num"), action="store", dest="burnin_num", type="integer", help="Number of burnin steps"),
|
|
|
8 make_option(c("--bychr"), action="store_true", dest="bychr", default=FALSE, help="Output chromosomes in separate files"),
|
|
|
9 make_option(c("--hp"), action="store_true", dest="hp", default=FALSE, help="Discourage state transition across chromosomes"),
|
|
|
10 make_option(c("--initial_states"), action="store", dest="initial_states", type="integer", default=NULL, help="Initial number of states"),
|
|
|
11 make_option(c("--log2"), action="store", dest="log2", type="double", default=NULL, help="log2 transformation"),
|
|
|
12 make_option(c("--maxerr"), action="store", dest="maxerr", type="double", default=NULL, help="Maximum standard deviation for the emission Gaussian distribution"),
|
|
|
13 make_option(c("--max_cell_type_clusters"), action="store", dest="max_cell_type_clusters", type="integer", default=NULL, help="Maximum number of cell type clusters allowed"),
|
|
|
14 make_option(c("--max_position_classes"), action="store", dest="max_position_classes", type="integer", default=NULL, help="Maximum number of position classes to be inferred"),
|
|
|
15 make_option(c("--max_states"), action="store", dest="max_states", type="double", default=NULL, help="Maximum number of states to be inferred"),
|
|
|
16 make_option(c("--mcmc_num"), action="store", dest="mcmc_num", type="integer", help="Number of maximization steps"),
|
|
|
17 make_option(c("--minerr"), action="store", dest="minerr", type="double", default=NULL, help="Minimum standard deviation for the emission Gaussian distribution"),
|
|
|
18 make_option(c("--norm"), action="store_true", dest="norm", default=FALSE, help="Standardize all datasets"),
|
|
|
19 make_option(c("--output_log"), action="store", dest="output_log", default=NULL, help="Output log file path"),
|
|
|
20 make_option(c("--prep_output_config"), action="store", dest="prep_output_config", help="prepMat output config file"),
|
|
|
21 make_option(c("--prior_concentration"), action="store", dest="prior_concentration", type="double", default=NULL, help="Prior concentration"),
|
|
|
22 make_option(c("--project_name"), action="store", dest="project_name", help="Outputs will have this base name"),
|
|
|
23 make_option(c("--rseed"), action="store", dest="rseed", type="integer", help="Seed for IDEAS model initialization"),
|
|
|
24 make_option(c("--save_ideas_log"), action="store", dest="save_ideas_log", default=NULL, help="Flag to save IDEAS process log"),
|
|
|
25 make_option(c("--script_dir"), action="store", dest="script_dir", help="R script source directory"),
|
|
|
26 make_option(c("--thread"), action="store", dest="thread", type="integer", help="Process threads"),
|
|
|
27 make_option(c("--tmp_dir"), action="store", dest="tmp_dir", help="Directory of bed files"),
|
|
|
28 make_option(c("--training_iterations"), action="store", dest="training_iterations", type="integer", default=NULL, help="Number of training iterations"),
|
|
|
29 make_option(c("--training_windows"), action="store", dest="training_windows", type="integer", default=NULL, help="Number of training iterations"),
|
|
|
30 make_option(c("--windows_bed"), action="store", dest="windows_bed", default=NULL, help="Bed file containing bed windows"),
|
|
166
|
31 make_option(c("--windows_config"), action="store", dest="windows_config", default=NULL, help="Windows positions by chroms config")
|
|
150
|
32 )
|
|
|
33
|
|
|
34 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
|
|
|
35 args <- parse_args(parser, positional_arguments=TRUE)
|
|
|
36 opt <- args$options
|
|
|
37
|
|
|
38 add_output_redirect <- function(cmd, save_ideas_log, output_log, default_log_name) {
|
|
|
39 if (is.null(save_ideas_log)) {
|
|
166
|
40 new_cmd = c(cmd, "&>>", default_log_name);
|
|
150
|
41 }else {
|
|
166
|
42 new_cmd = c(cmd, "&>>", output_log);
|
|
150
|
43 }
|
|
166
|
44 return(paste(new_cmd, collapse=" "));
|
|
150
|
45 }
|
|
|
46
|
|
|
47 combine_state <- function(parafiles, method="ward.D", mycut=0.9, pcut=1.0) {
|
|
|
48 X = NULL;
|
|
|
49 K = NULL;
|
|
|
50 I = NULL;
|
|
|
51 myheader = NULL;
|
|
|
52 p = NULL;
|
|
|
53 for(i in 1:length(parafiles)) {
|
|
|
54 x = fread(parafiles[i]);
|
|
|
55 t = max(which(is.na(x[1,])==F));
|
|
|
56 x = as.matrix(x[,1:t]);
|
|
|
57 if(i==1) {
|
|
|
58 myheader = colnames(x);
|
|
|
59 p = sqrt(9/4-2*(1-length(myheader))) - 3 / 2;
|
|
|
60 }
|
|
|
61 m = match(myheader[1:p+1], colnames(x)[1:p+1]);
|
|
|
62 v = NULL;
|
|
|
63 for(ii in 1:p) {
|
|
|
64 for(jj in 1:ii) {
|
|
|
65 a = max(m[ii],m[jj]);
|
|
|
66 b = min(m[ii],m[jj]);
|
|
|
67 v = c(v, a*(a+1)/2+b-a);
|
|
|
68 }
|
|
|
69 }
|
|
|
70 X = rbind(X, array(as.matrix(x[, c(1, 1+m, 1+p+v)]), dim=c(length(x) / (1+p+length(v)), 1 + p + length(v))));
|
|
|
71 K = c(K, dim(x)[1]);
|
|
|
72 I = c(I, rep(i, dim(x)[1]));
|
|
|
73 }
|
|
|
74 N = length(parafiles);
|
|
|
75 p = (sqrt(1 + dim(X)[2] * 8) - 3) / 2;
|
|
|
76 omycut = mycut;
|
|
|
77 mycut = round(length(parafiles) * mycut);
|
|
|
78 M = array(X[,1:p+1] / X[,1], dim=c(dim(X)[1], p));
|
|
|
79 V = array(0, dim=c(dim(X)[1] * p, p));
|
|
|
80 for(i in 1:dim(X)[1]) {
|
|
|
81 t = (i - 1) * p;
|
|
|
82 l = 1;
|
|
|
83 for(j in 1:p) {
|
|
|
84 for(k in 1:j) {
|
|
|
85 V[t+j, k] = V[t+k, j] = X[i,1+p+l] / X[i,1] - M[i,j] * M[i,k];
|
|
|
86 l = l + 1;
|
|
|
87 }
|
|
|
88 }
|
|
|
89 V[t+1:p,] = t(solve(chol(V[t+1:p,] + diag(1e-1,p))));
|
|
|
90 }
|
|
|
91 D = array(0, dim=rep(dim(X)[1],2));
|
|
|
92 for(i in 2:dim(X)[1]) {
|
|
|
93 for(j in 1:(i-1)) {
|
|
|
94 D[i,j] = D[j,i] = sqrt((sum((V[(i-1)*p+1:p,]%*%(M[i,]-M[j,]))^2) + sum((V[(j-1)*p+1:p,]%*%(M[i,]-M[j,]))^2)));
|
|
|
95 }
|
|
|
96 }
|
|
|
97 MM = NULL;
|
|
|
98 kk = NULL;
|
|
|
99 for(i in 1:N) {
|
|
|
100 t = 1:K[i];
|
|
|
101 if(i > 1) {
|
|
|
102 t = t + sum(K[1:(i-1)]);
|
|
|
103 }
|
|
|
104 t = (1:dim(D)[1])[-t];
|
|
|
105 h = hclust(as.dist(D[t,t]), method=method);
|
|
|
106 k = -1;
|
|
|
107 tM = NULL;
|
|
|
108 for(j in min(K):(min(length(t), max(K)*2))) {
|
|
|
109 m = cutree(h,k=j);
|
|
|
110 tt = NULL;
|
|
|
111 for(l in 1:j) {
|
|
|
112 tt[l] = length(unique(I[t[which(m==l)]]));
|
|
|
113 }
|
|
|
114 tk = length(which(tt>=mycut));
|
|
|
115 if(tk > k) {
|
|
|
116 k = tk;
|
|
|
117 tM = make_parameter(1:j, I[t], m, mycut, X[t,]);
|
|
|
118 } else if(tk < k) {
|
|
|
119 break;
|
|
|
120 }
|
|
|
121 }
|
|
|
122 kk[i] = k;
|
|
|
123 MM = rbind(MM, cbind(i, tM));
|
|
|
124 }
|
|
|
125 mysel = median(kk);
|
|
|
126 h = hclust(as.dist(D), method=method);
|
|
|
127 rt = rep(0, max(K)*2);
|
|
|
128 k = -1;
|
|
|
129 for(i in min(K):min(dim(D)[1], max(K)*2)) {
|
|
|
130 m = cutree(h,k=i);
|
|
|
131 tt = NULL;
|
|
|
132 for(j in 1:i) {
|
|
|
133 tt[j] = length(unique(I[which(m==j)]));
|
|
|
134 }
|
|
|
135 tk = length(which(tt>=mycut));
|
|
|
136 if(tk==mysel | tk<k) {
|
|
|
137 break;
|
|
|
138 }
|
|
|
139 k = tk;
|
|
|
140 rt[i] = length(which(tt>=mycut));
|
|
|
141 }
|
|
|
142 mysel = max(k,tk);
|
|
|
143 m = cutree(h, k=mysel);
|
|
|
144 nn = NULL;
|
|
|
145 for(i in 1:mysel) {
|
|
|
146 t = which(m==i);
|
|
|
147 nn[i] = sum(X[t,1]);
|
|
|
148 }
|
|
|
149 oo = order(nn, decreasing=T);
|
|
|
150 rt = make_parameter(oo, I, m, mycut, X);
|
|
|
151 onstate = max(rt[,1]) + 1;
|
|
|
152 ooo = NULL;
|
|
|
153 for(i in oo) {
|
|
|
154 t = which(m==i);
|
|
|
155 if(length(unique(I[t])) >= mycut) {
|
|
|
156 ooo = c(ooo, i);
|
|
|
157 }
|
|
|
158 }
|
|
|
159 d = NULL;
|
|
|
160 for(i in 1:N) {
|
|
|
161 d = rbind(d, compare_two(rt, MM[MM[,1]==i,-1])[1:onstate]);
|
|
|
162 }
|
|
|
163 dd = array(cutree(hclust(dist(c(d))), k=2), dim=dim(d));
|
|
|
164 kk = table(c(dd));
|
|
|
165 kk = which(as.integer(kk)==max(as.integer(kk)))[1];
|
|
|
166 pp = apply(dd, 2, function(x){length(which(x!=kk))/length(x)});
|
|
|
167 pp0 = apply(d, 2, function(x){length(which(x>0.5))/length(x)});
|
|
|
168 pp[pp0<pp] = pp0[pp0<pp];
|
|
|
169 t = which(pp > pcut);
|
|
|
170 if(length(t) > 0) {
|
|
|
171 j = 0;
|
|
|
172 nrt = NULL;
|
|
|
173 for(i in (1:onstate-1)[-t]) {
|
|
|
174 nrt = rbind(nrt, cbind(j, rt[rt[,1]==i,-1]));
|
|
|
175 j = j + 1;
|
|
|
176 }
|
|
|
177 rt = nrt;
|
|
|
178 ooo = ooo[-t];
|
|
|
179 }
|
|
|
180 nrt = NULL;
|
|
|
181 for(i in 0:max(rt[,1])) {
|
|
|
182 t = which(rt[,1]==i);
|
|
|
183 nrt = rbind(nrt, apply(array(rt[t,], dim=c(length(t), dim(rt)[2])), 2, sum)[-1]);
|
|
|
184 }
|
|
|
185 rt = nrt;
|
|
|
186 colnames(rt) = myheader;
|
|
|
187 O = NULL;
|
|
|
188 Ip = NULL;
|
|
|
189 Xp = NULL;
|
|
|
190 k = 0;
|
|
|
191 for(i in 1:length(parafiles)) {
|
|
|
192 str = gsub(".para", ".profile", parafiles[i]);
|
|
|
193 p = as.matrix(read.table(str));
|
|
|
194 u = array(0, dim=c(dim(p)[1], length(ooo)));
|
|
|
195 for(j in 1:length(ooo)) {
|
|
|
196 t = which(m[k+1:K[i]] == ooo[j]);
|
|
|
197 u[,j] = apply(array(p[,1+t], dim=c(dim(p)[1], length(t))), 1, sum);
|
|
|
198 }
|
|
|
199 k = k + K[i];
|
|
|
200 u = u / (apply(u, 1, sum) + 1e-10);
|
|
|
201 Xp = rbind(Xp, cbind(p[,1], u));
|
|
|
202 Ip = c(Ip, rep(i,dim(u)[1]));
|
|
|
203 }
|
|
|
204 hp = hclust(dist(((Xp[,-1]+min(1e-3, min(Xp[,-1][Xp[,-1]>0]))))), method=method);
|
|
|
205 ocut = min(mycut/2, length(parafiles)/2);
|
|
|
206 t = range(as.integer(table(Ip)));
|
|
|
207 Kp = NULL;
|
|
|
208 for(i in t[1]:(t[2]*2)) {
|
|
|
209 m = cutree(hp, k=i);
|
|
|
210 tt = table(Ip,m);
|
|
|
211 ll = apply(tt, 2, function(x){length(which(x>0))});
|
|
|
212 Kp = c(Kp, length(which(ll>=ocut)));
|
|
|
213 }
|
|
|
214 oN = (t[1]:(t[2]*2))[which(Kp==max(Kp))[1]];
|
|
|
215 m = cutree(hp, k=oN);
|
|
|
216 tt = table(Ip,m);
|
|
|
217 ll = apply(tt, 2, function(x){length(which(x>0))});
|
|
|
218 tt = which(ll>=ocut);
|
|
|
219 for(i in tt) {
|
|
|
220 t = which(m==i);
|
|
|
221 O = rbind(O, c(sum(Xp[t, 1]), apply(array(Xp[t,-1]*Xp[t,1], dim=c(length(t), dim(Xp)[2]-1)), 2, sum)/sum(Xp[t, 1])));
|
|
|
222 }
|
|
|
223 nrt = NULL;
|
|
|
224 nrt$para = rt;
|
|
|
225 nrt$profile = O;
|
|
|
226 return(nrt);
|
|
|
227 }
|
|
|
228
|
|
|
229 compare_two <- function(n, m) {
|
|
|
230 NN = get_mean(n);
|
|
|
231 MM = get_mean(m);
|
|
|
232 p = (-3 + sqrt(9 + 8 * (dim(n)[2] - 2))) / 2;
|
|
|
233 dd = NULL;
|
|
|
234 for (i in 1:dim(NN)[1]) {
|
|
|
235 dd[i] = min(apply(array(MM[,1:p], dim=c(dim(MM)[1],p)), 1, function(x){sqrt(sum((x-NN[i,1:p])^2))}));
|
|
|
236 }
|
|
|
237 for (i in 1:dim(MM)[1]) {
|
|
|
238 dd[i+dim(NN)[1]] = min(apply(array(NN[,1:p], dim=c(dim(NN)[1],p)), 1, function(x){sqrt(sum((x-MM[i,1:p])^2))}));
|
|
|
239 }
|
|
|
240 return(dd);
|
|
|
241 }
|
|
|
242
|
|
166
|
243 get_base_cmd <- function(prep_output_config, windows_bed, training_iterations, bychr, hp, norm, log2,
|
|
|
244 max_states, initial_states, max_position_classes, max_cell_type_clusters, prior_concentration,
|
|
|
245 burnin_num, mcmc_num, minerr, maxerr, rseed, thread) {
|
|
|
246 base_cmd = paste("ideas", prep_output_config, sep=" ");
|
|
|
247 if (!is.null(windows_bed)) {
|
|
|
248 base_cmd = paste(base_cmd, windows_bed, sep=" ");
|
|
|
249 }
|
|
|
250 if (!is.null(training_iterations)) {
|
|
|
251 base_cmd = paste(base_cmd, "-impute none", sep=" ");
|
|
|
252 }
|
|
|
253 if (bychr) {
|
|
|
254 base_cmd = paste(base_cmd, "-bychr", sep=" ");
|
|
|
255 }
|
|
|
256 if (hp) {
|
|
|
257 base_cmd = paste(base_cmd, "-hp", sep=" ");
|
|
|
258 }
|
|
|
259 if (norm) {
|
|
|
260 base_cmd = paste(base_cmd, "-norm", sep=" ");
|
|
|
261 }
|
|
|
262 if (!is.null(log2)) {
|
|
|
263 base_cmd = paste(base_cmd, "-log2", log2, sep=" ");
|
|
|
264 }
|
|
|
265 if (!is.null(max_states)) {
|
|
|
266 base_cmd = paste(base_cmd, "-G", max_states, sep=" ");
|
|
|
267 }
|
|
|
268 if (!is.null(initial_states)) {
|
|
|
269 base_cmd = paste(base_cmd, "-C", initial_states, sep=" ");
|
|
|
270 }
|
|
|
271 if (!is.null(max_position_classes)) {
|
|
|
272 base_cmd = paste(base_cmd, "-P", max_position_classes, sep=" ");
|
|
|
273 }
|
|
|
274 if (!is.null(max_cell_type_clusters)) {
|
|
|
275 base_cmd = paste(base_cmd, "-K", max_cell_type_clusters, sep=" ");
|
|
|
276 }
|
|
|
277 if (!is.null(prior_concentration)) {
|
|
|
278 base_cmd = paste(base_cmd, "-A", prior_concentration, sep=" ");
|
|
|
279 }
|
|
|
280 base_cmd = paste(base_cmd, "-sample", burnin_num, mcmc_num, sep=" ");
|
|
|
281 if (!is.null(minerr)) {
|
|
|
282 base_cmd = paste(base_cmd, "-minerr", minerr, sep=" ");
|
|
|
283 }
|
|
|
284 if (!is.null(maxerr)) {
|
|
|
285 base_cmd = paste(base_cmd, "-maxerr", maxerr, sep=" ");
|
|
|
286 }
|
|
|
287 base_cmd = paste(base_cmd, "-rseed", rseed, sep=" ");
|
|
|
288 base_cmd = paste(base_cmd, "-thread", thread, sep=" ");
|
|
|
289 return(base_cmd);
|
|
|
290 }
|
|
|
291
|
|
150
|
292 get_mean <- function(n) {
|
|
|
293 N = NULL;
|
|
|
294 for(i in sort(unique(n[,1]))) {
|
|
|
295 t = which(n[,1]==i);
|
|
|
296 N = rbind(N, apply(array(n[t,], dim=c(length(t), dim(n)[2])), 2, sum)[-1]);
|
|
|
297 }
|
|
|
298 NN = N[,-1] / N[,1];
|
|
|
299 return(array(NN, dim=c(length(NN)/(dim(n)[2]-2), dim(n)[2]-2)));
|
|
|
300 }
|
|
|
301
|
|
166
|
302 get_post_training_base_cmd <- function(base_cmd, para) {
|
|
|
303 # Change base_cmd due to training mode.
|
|
|
304 base_cmd_items = as.list(strsplit(base_cmd[1], split=" ", fixed=TRUE))[[1]];
|
|
|
305 if (length(which(base_cmd_items == "-G")) == 0) {
|
|
|
306 base_cmd_items = c(base_cmd_items, "-G", length(para)-1);
|
|
|
307 } else {
|
|
|
308 tt = which(base_cmd_items == "-G");
|
|
|
309 base_cmd_items[tt + 1] = length(para)-1;
|
|
|
310 }
|
|
|
311 tt = which(base_cmd_items == '-C');
|
|
|
312 if(length(tt) > 0) {
|
|
|
313 base_cmd_items = base_cmd_items[-c(tt, tt+1)];
|
|
|
314 }
|
|
|
315 base_cmd = paste(base_cmd_items, collapse=" ");
|
|
|
316 return(base_cmd);
|
|
|
317 }
|
|
|
318
|
|
|
319 get_windows_by_chrom <- function(windows_config) {
|
|
|
320 if (is.null(windows_config)) {
|
|
|
321 windows_by_chrom = NULL;
|
|
|
322 } else {
|
|
|
323 fh = file(windows_config, "r");
|
|
|
324 windows_by_chrom = readLines(fh);
|
|
|
325 close(fh);
|
|
|
326 }
|
|
|
327 return(windows_by_chrom);
|
|
|
328 }
|
|
|
329
|
|
150
|
330 make_parameter <- function(myorder, id, mem, mycut, para) {
|
|
|
331 rt = NULL;
|
|
|
332 j = 0;
|
|
|
333 for(i in myorder) {
|
|
|
334 t = which(mem==i);
|
|
|
335 if (length(unique(id[t])) >= mycut) {
|
|
|
336 rt = rbind(rt, cbind(j, array(para[t,], dim=c(length(t), dim(para)[2]))));
|
|
|
337 j = j + 1;
|
|
|
338 }
|
|
|
339 }
|
|
|
340 return(rt);
|
|
|
341 }
|
|
|
342
|
|
166
|
343 remove_files <- function(path, pattern) {
|
|
|
344 files = list.files(path=path, pattern=pattern);
|
|
|
345 for (f in files) {
|
|
|
346 unlink(f);
|
|
|
347 }
|
|
|
348 }
|
|
|
349
|
|
150
|
350 run_cmd <- function(cmd, save_ideas_log, output_log, default_log_name) {
|
|
|
351 rc = system(cmd);
|
|
|
352 if (rc != 0) {
|
|
|
353 if (is.null(save_ideas_log)) {
|
|
|
354 file.rename(default_log_name, output_log);
|
|
|
355 }
|
|
|
356 quit(rc);
|
|
|
357 }
|
|
|
358 }
|
|
|
359
|
|
|
360 default_log_name = "ideas_log.txt";
|
|
166
|
361 windows_by_chrom = get_windows_by_chrom(opt$windows_config);
|
|
|
362 base_cmd = get_base_cmd(opt$prep_output_config, opt$windows_bed, opt$training_iterations, opt$bychr,
|
|
|
363 opt$hp, opt$norm, opt$log2, opt$max_states, opt$initial_states, opt$max_position_classes,
|
|
|
364 opt$max_cell_type_clusters, opt$prior_concentration, opt$burnin_num, opt$mcmc_num, opt$minerr,
|
|
|
365 opt$maxerr, opt$rseed, opt$thread);
|
|
150
|
366 output_base_name = opt$project_name;
|
|
|
367
|
|
|
368 if (is.null(opt$training_iterations)) {
|
|
166
|
369 # Not performing training.
|
|
|
370 if (is.null(windows_by_chrom)) {
|
|
|
371 # Not performing windows by chromosome.
|
|
|
372 output_name = output_base_name;
|
|
|
373 cmd = paste(base_cmd, "-o", output_name, sep=" ");
|
|
|
374 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
375 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
376 } else {
|
|
|
377 # Performing windows by chromosome.
|
|
|
378 for (i in 1:length(windows_by_chrom)) {
|
|
|
379 line = windows_by_chrom[i];
|
|
|
380 items = strsplit(line, " ")[[1]];
|
|
|
381 chrom = items[1];
|
|
|
382 window_start = items[2];
|
|
|
383 window_end = items[3];
|
|
|
384 output_name = paste(output_base_name,, chrom, sep=".");
|
|
|
385 cmd = paste(base_cmd, "-inv", window_start, window_end, sep=" ");
|
|
|
386 cmd = paste(cmd, "-o", output_name, sep=" ");
|
|
|
387 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
388 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
389 }
|
|
|
390 }
|
|
150
|
391 } else {
|
|
166
|
392 # performing training.
|
|
|
393 output_para0 = paste(output_base_name, "para0", sep=".");
|
|
|
394 output_profile0 = paste(output_base_name, "profile0", sep=".");
|
|
150
|
395 for (i in 1:opt$training_iterations) {
|
|
166
|
396 cmd = paste(base_cmd, "-o", paste(output_base_name, ".tmp.", i, sep=""), sep=" ");
|
|
|
397 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
398 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
150
|
399 }
|
|
166
|
400 tpara = combine_state(paste(output_base_name, "tmp", (1:opt$training_iterations), "para", sep="."), mycut=0.5);
|
|
|
401 write.table(tpara$profile, output_profile0, quote=F, row.names=F, col.names=F);
|
|
150
|
402 para = tpara$para;
|
|
|
403 para = apply(para, 1, function(x){paste(x, collapse=" ")});
|
|
166
|
404 para = c(readLines(paste(output_base_name, "tmp", "1", "para", sep="."), n=1), para);
|
|
150
|
405 writeLines(para, output_para0);
|
|
166
|
406 # Now run IDEAS based on the files produced during training.
|
|
|
407 base_cmd = get_post_training_base_cmd(base_cmd, para);
|
|
|
408 base_cmd = paste(base_cmd, "-otherpara", output_para0[[1]], output_profile0[[1]], sep=" ");
|
|
|
409 if (is.null(windows_by_chrom)) {
|
|
|
410 cmd = c(base_cmd, "-o", output_base_name);
|
|
|
411 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
412 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
152
|
413 } else {
|
|
166
|
414 # Performing windows by chromosome.
|
|
|
415 if (length(windows_by_chrom) == 1) {
|
|
|
416 output_name = paste(output_base_name, i, sep=".");
|
|
|
417 cmd = c(base_cmd, "-o", output_name);
|
|
|
418 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
419 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
420 } else {
|
|
|
421 for (i in 1:length(windows_by_chrom)) {
|
|
|
422 line = windows_by_chrom[i];
|
|
|
423 items = strsplit(line, " ")[[1]];
|
|
|
424 chrom = items[[1]];
|
|
|
425 window_start = items[[2]];
|
|
|
426 window_end = items[[3]];
|
|
|
427 cmd = paste(base_cmd, "-inv", window_start, window_end, sep=" ");
|
|
|
428 output_name = paste(output_base_name, chrom, sep=".");
|
|
|
429 cmd = paste(cmd, "-o", output_name, sep=" ");
|
|
|
430 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
431 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
|
|
|
432 }
|
|
|
433 }
|
|
152
|
434 }
|
|
166
|
435 # Remove temporary outputs.
|
|
|
436 remove_files(path=".", pattern="tmp");
|
|
150
|
437 }
|
