Mercurial > repos > greg > ideas

comparison ideas.R @ 150:3762c27d820a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 12 Jan 2018 11:24:43 -0500
parents
children 26c26cb32137
comparison
equal deleted inserted replaced
149:a80b76535243 150:3762c27d820a
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("--output_txt_dir"), action="store", dest="output_txt_dir", help="Directory for output txt files"),
21 make_option(c("--prep_output_config"), action="store", dest="prep_output_config", help="prepMat output config file"),
22 make_option(c("--prior_concentration"), action="store", dest="prior_concentration", type="double", default=NULL, help="Prior concentration"),
23 make_option(c("--project_name"), action="store", dest="project_name", help="Outputs will have this base name"),
24 make_option(c("--rseed"), action="store", dest="rseed", type="integer", help="Seed for IDEAS model initialization"),
25 make_option(c("--save_ideas_log"), action="store", dest="save_ideas_log", default=NULL, help="Flag to save IDEAS process log"),
26 make_option(c("--script_dir"), action="store", dest="script_dir", help="R script source directory"),
27 make_option(c("--thread"), action="store", dest="thread", type="integer", help="Process threads"),
28 make_option(c("--tmp_dir"), action="store", dest="tmp_dir", help="Directory of bed files"),
29 make_option(c("--training_iterations"), action="store", dest="training_iterations", type="integer", default=NULL, help="Number of training iterations"),
30 make_option(c("--training_windows"), action="store", dest="training_windows", type="integer", default=NULL, help="Number of training iterations"),
31 make_option(c("--windows_bed"), action="store", dest="windows_bed", default=NULL, help="Bed file containing bed windows"),
32 make_option(c("--window_end"), action="store", dest="window_end", type="integer", default=NULL, help="Windows positions by chromosome end value"),
33 make_option(c("--window_start"), action="store", dest="window_start", type="integer", default=NULL, help="Windows positions by chromosome start value")
34 )
35
36 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
37 args <- parse_args(parser, positional_arguments=TRUE)
38 opt <- args$options
39
40 add_output_redirect <- function(cmd, save_ideas_log, output_log, default_log_name) {
41 if (is.null(save_ideas_log)) {
42 cmd = paste(cmd, "&>>", default_log_name, sep=" ");
43 }else {
44 cmd = paste(cmd, "&>>", output_log, sep=" ");
45 }
46 return(cmd);
47 }
48
49 combine_state <- function(parafiles, method="ward.D", mycut=0.9, pcut=1.0) {
50 X = NULL;
51 K = NULL;
52 I = NULL;
53 myheader = NULL;
54 p = NULL;
55 for(i in 1:length(parafiles)) {
56 x = fread(parafiles[i]);
57 t = max(which(is.na(x[1,])==F));
58 x = as.matrix(x[,1:t]);
59 if(i==1) {
60 myheader = colnames(x);
61 p = sqrt(9/4-2*(1-length(myheader))) - 3 / 2;
62 }
63 m = match(myheader[1:p+1], colnames(x)[1:p+1]);
64 v = NULL;
65 for(ii in 1:p) {
66 for(jj in 1:ii) {
67 a = max(m[ii],m[jj]);
68 b = min(m[ii],m[jj]);
69 v = c(v, a*(a+1)/2+b-a);
70 }
71 }
72 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))));
73 K = c(K, dim(x)[1]);
74 I = c(I, rep(i, dim(x)[1]));
75 }
76 N = length(parafiles);
77 p = (sqrt(1 + dim(X)[2] * 8) - 3) / 2;
78 omycut = mycut;
79 mycut = round(length(parafiles) * mycut);
80 M = array(X[,1:p+1] / X[,1], dim=c(dim(X)[1], p));
81 V = array(0, dim=c(dim(X)[1] * p, p));
82 for(i in 1:dim(X)[1]) {
83 t = (i - 1) * p;
84 l = 1;
85 for(j in 1:p) {
86 for(k in 1:j) {
87 V[t+j, k] = V[t+k, j] = X[i,1+p+l] / X[i,1] - M[i,j] * M[i,k];
88 l = l + 1;
89 }
90 }
91 V[t+1:p,] = t(solve(chol(V[t+1:p,] + diag(1e-1,p))));
92 }
93 D = array(0, dim=rep(dim(X)[1],2));
94 for(i in 2:dim(X)[1]) {
95 for(j in 1:(i-1)) {
96 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)));
97 }
98 }
99 MM = NULL;
100 kk = NULL;
101 for(i in 1:N) {
102 t = 1:K[i];
103 if(i > 1) {
104 t = t + sum(K[1:(i-1)]);
105 }
106 t = (1:dim(D)[1])[-t];
107 h = hclust(as.dist(D[t,t]), method=method);
108 k = -1;
109 tM = NULL;
110 for(j in min(K):(min(length(t), max(K)*2))) {
111 m = cutree(h,k=j);
112 tt = NULL;
113 for(l in 1:j) {
114 tt[l] = length(unique(I[t[which(m==l)]]));
115 }
116 tk = length(which(tt>=mycut));
117 if(tk > k) {
118 k = tk;
119 tM = make_parameter(1:j, I[t], m, mycut, X[t,]);
120 } else if(tk < k) {
121 break;
122 }
123 }
124 kk[i] = k;
125 MM = rbind(MM, cbind(i, tM));
126 }
127 mysel = median(kk);
128 h = hclust(as.dist(D), method=method);
129 rt = rep(0, max(K)*2);
130 k = -1;
131 for(i in min(K):min(dim(D)[1], max(K)*2)) {
132 m = cutree(h,k=i);
133 tt = NULL;
134 for(j in 1:i) {
135 tt[j] = length(unique(I[which(m==j)]));
136 }
137 tk = length(which(tt>=mycut));
138 if(tk==mysel | tk<k) {
139 break;
140 }
141 k = tk;
142 rt[i] = length(which(tt>=mycut));
143 }
144 mysel = max(k,tk);
145 m = cutree(h, k=mysel);
146 nn = NULL;
147 for(i in 1:mysel) {
148 t = which(m==i);
149 nn[i] = sum(X[t,1]);
150 }
151 oo = order(nn, decreasing=T);
152 rt = make_parameter(oo, I, m, mycut, X);
153 onstate = max(rt[,1]) + 1;
154 ooo = NULL;
155 for(i in oo) {
156 t = which(m==i);
157 if(length(unique(I[t])) >= mycut) {
158 ooo = c(ooo, i);
159 }
160 }
161 d = NULL;
162 for(i in 1:N) {
163 d = rbind(d, compare_two(rt, MM[MM[,1]==i,-1])[1:onstate]);
164 }
165 dd = array(cutree(hclust(dist(c(d))), k=2), dim=dim(d));
166 kk = table(c(dd));
167 kk = which(as.integer(kk)==max(as.integer(kk)))[1];
168 pp = apply(dd, 2, function(x){length(which(x!=kk))/length(x)});
169 pp0 = apply(d, 2, function(x){length(which(x>0.5))/length(x)});
170 pp[pp0<pp] = pp0[pp0<pp];
171 t = which(pp > pcut);
172 if(length(t) > 0) {
173 j = 0;
174 nrt = NULL;
175 for(i in (1:onstate-1)[-t]) {
176 nrt = rbind(nrt, cbind(j, rt[rt[,1]==i,-1]));
177 j = j + 1;
178 }
179 rt = nrt;
180 ooo = ooo[-t];
181 }
182 nrt = NULL;
183 for(i in 0:max(rt[,1])) {
184 t = which(rt[,1]==i);
185 nrt = rbind(nrt, apply(array(rt[t,], dim=c(length(t), dim(rt)[2])), 2, sum)[-1]);
186 }
187 rt = nrt;
188 colnames(rt) = myheader;
189 O = NULL;
190 Ip = NULL;
191 Xp = NULL;
192 k = 0;
193 for(i in 1:length(parafiles)) {
194 str = gsub(".para", ".profile", parafiles[i]);
195 p = as.matrix(read.table(str));
196 u = array(0, dim=c(dim(p)[1], length(ooo)));
197 for(j in 1:length(ooo)) {
198 t = which(m[k+1:K[i]] == ooo[j]);
199 u[,j] = apply(array(p[,1+t], dim=c(dim(p)[1], length(t))), 1, sum);
200 }
201 k = k + K[i];
202 u = u / (apply(u, 1, sum) + 1e-10);
203 Xp = rbind(Xp, cbind(p[,1], u));
204 Ip = c(Ip, rep(i,dim(u)[1]));
205 }
206 hp = hclust(dist(((Xp[,-1]+min(1e-3, min(Xp[,-1][Xp[,-1]>0]))))), method=method);
207 ocut = min(mycut/2, length(parafiles)/2);
208 t = range(as.integer(table(Ip)));
209 Kp = NULL;
210 for(i in t[1]:(t[2]*2)) {
211 m = cutree(hp, k=i);
212 tt = table(Ip,m);
213 ll = apply(tt, 2, function(x){length(which(x>0))});
214 Kp = c(Kp, length(which(ll>=ocut)));
215 }
216 oN = (t[1]:(t[2]*2))[which(Kp==max(Kp))[1]];
217 m = cutree(hp, k=oN);
218 tt = table(Ip,m);
219 ll = apply(tt, 2, function(x){length(which(x>0))});
220 tt = which(ll>=ocut);
221 for(i in tt) {
222 t = which(m==i);
223 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])));
224 }
225 nrt = NULL;
226 nrt$para = rt;
227 nrt$profile = O;
228 return(nrt);
229 }
230
231 compare_two <- function(n, m) {
232 NN = get_mean(n);
233 MM = get_mean(m);
234 p = (-3 + sqrt(9 + 8 * (dim(n)[2] - 2))) / 2;
235 dd = NULL;
236 for (i in 1:dim(NN)[1]) {
237 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))}));
238 }
239 for (i in 1:dim(MM)[1]) {
240 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))}));
241 }
242 return(dd);
243 }
244
245 get_mean <- function(n) {
246 N = NULL;
247 for(i in sort(unique(n[,1]))) {
248 t = which(n[,1]==i);
249 N = rbind(N, apply(array(n[t,], dim=c(length(t), dim(n)[2])), 2, sum)[-1]);
250 }
251 NN = N[,-1] / N[,1];
252 return(array(NN, dim=c(length(NN)/(dim(n)[2]-2), dim(n)[2]-2)));
253 }
254
255 make_parameter <- function(myorder, id, mem, mycut, para) {
256 rt = NULL;
257 j = 0;
258 for(i in myorder) {
259 t = which(mem==i);
260 if (length(unique(id[t])) >= mycut) {
261 rt = rbind(rt, cbind(j, array(para[t,], dim=c(length(t), dim(para)[2]))));
262 j = j + 1;
263 }
264 }
265 return(rt);
266 }
267
268 run_cmd <- function(cmd, save_ideas_log, output_log, default_log_name) {
269 cat("\n\n >>>>> cmd:\n", cmd, "\n\n");
270 rc = system(cmd);
271 if (rc != 0) {
272 if (is.null(save_ideas_log)) {
273 file.rename(default_log_name, output_log);
274 }
275 quit(rc);
276 }
277 }
278
279 default_log_name = "ideas_log.txt";
280 output_base_name = opt$project_name;
281 cmd = paste("ideas", opt$prep_output_config, sep=" ");
282 if (!is.null(opt$windows_bed)) {
283 cmd = paste(cmd, opt$windows_bed, sep=" ");
284 }
285 if (!is.null(opt$training_iterations)) {
286 cmd = paste(cmd, "-impute none", sep=" ");
287 }
288 if (opt$bychr) {
289 cmd = paste(cmd, "-bychr", sep=" ");
290 }
291 if (opt$hp) {
292 cmd = paste(cmd, "-hp", sep=" ");
293 }
294 if (opt$norm) {
295 cmd = paste(cmd, "-norm", sep=" ");
296 }
297 if (!is.null(opt$window_start) && !is.null(opt$window_end)) {
298 cmd = paste(cmd, "-inv", opt$window_start, opt$window_end, sep=" ");
299 }
300 if (!is.null(opt$log2)) {
301 cmd = paste(cmd, "-log2", opt$log2, sep=" ");
302 }
303 if (!is.null(opt$max_states)) {
304 cmd = paste(cmd, "-G", opt$max_states, sep=" ");
305 }
306 if (!is.null(opt$initial_states)) {
307 cmd = paste(cmd, "-C", opt$initial_states, sep=" ");
308 }
309 if (!is.null(opt$max_position_classes)) {
310 cmd = paste(cmd, "-P", opt$max_position_classes, sep=" ");
311 }
312 if (!is.null(opt$max_cell_type_clusters)) {
313 cmd = paste(cmd, "-K", opt$max_cell_type_clusters, sep=" ");
314 }
315 if (!is.null(opt$prior_concentration)) {
316 cmd = paste(cmd, "-A", opt$prior_concentration, sep=" ");
317 }
318 cmd = paste(cmd, "-sample", opt$burnin_num, opt$mcmc_num, sep=" ");
319 if (!is.null(opt$minerr)) {
320 cmd = paste(cmd, "-minerr", opt$minerr, sep=" ");
321 }
322 if (!is.null(opt$maxerr)) {
323 cmd = paste(cmd, "-maxerr", opt$maxerr, sep=" ");
324 }
325 cmd = paste(cmd, "-rseed", opt$rseed, sep=" ");
326 cmd = paste(cmd, "-thread", opt$thread, sep=" ");
327
328 if (is.null(opt$training_iterations)) {
329 cmd = paste(cmd, "-o", output_base_name, sep=" ");
330 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
331 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
332 } else {
333 for (i in 1:opt$training_iterations) {
334 cmd = paste(cmd, "-o", paste(output_base_name, ".tmp.", i, sep=""), sep=" ");
335 cmd = add_output_redirect(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
336 run_cmd(cmd, opt$save_ideas_log, opt$output_log, default_log_name);
337 }
338 tpara = combine_state(paste(output_base_name, ".tmp.", (1:opt$training_iterations), ".para", sep=""), mycut=0.5);
339 para = tpara$para;
340 write.table(tpara$profile, paste(output_base_name, ".profile0", sep=""), quote=F, row.names=F, col.names=F);
341 para = apply(para, 1, function(x){paste(x, collapse=" ")});
342 para = c(readLines(paste(output_base_name, ".tmp.1.para", sep=""), n=1), para);
343 output_para0 = paste(output_base_name, ".para0", sep="");
344 writeLines(para, output_para0);
345 }