Mercurial > repos > greg > ideas

comparison create_heatmap.R @ 118:f8a7bf4ca1a7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 17 Nov 2017 11:37:33 -0500
parents c89bf5bf4661
children 755d4a3754d2
comparison
equal deleted inserted replaced
117:42481a73f73d 118:f8a7bf4ca1a7
9 9
10 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 10 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
11 args <- parse_args(parser, positional_arguments=TRUE) 11 args <- parse_args(parser, positional_arguments=TRUE)
12 opt <- args$options 12 opt <- args$options
13 13
14 create_heatmap<-function(data_matrix, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) { 14 state_color <- function(statemean, markcolor=NULL)
15 k = dim(data_matrix)[2]; 15 {
16 l = dim(data_matrix)[1]; 16 if(length(markcolor) == 0) {
17 markcolor = rep("", dim(statemean)[2]);
18 markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:dim(statemean)[2]-1) / dim(statemean)[2], 1, 1);
19 markcolor = t(col2rgb(markcolor));
20 }
21
22 rg = apply(statemean, 1, range);
23 mm = NULL;
24 for(i in 1:dim(statemean)[1]) {
25 mm = rbind(mm, (statemean[i,] - rg[1,i] + 1e-10) / (rg[2,i] - rg[1,i] + 1e-10));
26 }
27 mm = mm^6;
28 if(dim(mm)[2] > 1) {
29 mm = mm / (apply(mm, 1, sum) + 1e-10);
30 }
31 mycol = mm%*%markcolor;
32 s = apply(statemean, 1, max);
33 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
34 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
35 h[,2] = h[,2] * s;
36 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
37 rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}) ,h);
38 return(rt);
39 }
40
41 create_heatmap<-function(data_frame, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) {
42 k = dim(data_frame)[2];
43 l = dim(data_frame)[1];
17 p = (sqrt(9 + 8 * (k - 1)) - 3) / 2; 44 p = (sqrt(9 + 8 * (k - 1)) - 3) / 2;
18 m = as.matrix(data_matrix[,1+1:p]/data_matrix[,1]); 45 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
19 colnames(m) = colnames(data_matrix)[1+1:p]; 46 colnames(data_matrix) = colnames(data_frame)[1+1:p];
20 marks = colnames(m); 47 marks = colnames(data_matrix);
21 rownames(m) = paste(1:l-1," (", round(data_matrix[,1] / sum(data_matrix[,1]) * 10000) / 100, "%)", sep=""); 48 rownames(data_matrix) = paste(1:l," (", round(data_frame[,1] / sum(data_frame[,1]) * 10000) / 100, "%)", sep="");
22 pdf(file=opt$output); 49 pdf(file=opt$output);
23 par(mar=c(6,1,1,6)); 50 par(mar=c(6,1,1,6));
24 rg = range(m); 51 rg = range(data_matrix);
25 colors = 0:100/100*(rg[2]-rg[1])+rg[1]; 52 colors = 0:100/100*(rg[2]-rg[1])+rg[1];
26 my_palette = colorRampPalette(cols)(n=100); 53 my_palette = colorRampPalette(cols)(n=100);
27 defpalette = palette(my_palette); 54 defpalette = palette(my_palette);
28 55
29 plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F); 56 plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
30 axis(1, at=1:p-0.5, labels=colnames(m), las=2); 57 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
31 axis(4, at=1:l-0.5, labels=rownames(m), las=2); 58 axis(4, at=1:l-0.5, labels=rownames(data_matrix), las=2);
32 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)); 59 rect(rep(1:p-1,l), rep(1:l-1,each=p), rep(1:p,l), rep(1:l,each=p), col=round((t(data_matrix)-rg[1])/(rg[2]-rg[1])*100));
33 60
34 if(length(statecolor)==0) { 61 if(length(statecolor)==0) {
35 if(length(markcolor)==0) { 62 if(length(markcolor)==0) {
36 markcolor=t(col2rgb(terrain.colors(ceiling(p))[1:p])); 63 markcolor=t(col2rgb(terrain.colors(ceiling(p))[1:p]));
37 for(i in 1:length(marks)) { 64 for(i in 1:length(marks)) {
74 if(regexpr("ctcf",tolower(marks[i]))>0) { 101 if(regexpr("ctcf",tolower(marks[i]))>0) {
75 markcolor[i,]=c(200,0,250); 102 markcolor[i,]=c(200,0,250);
76 } 103 }
77 } 104 }
78 } 105 }
79 statecolor = stateColor(m, markcolor)[,2]; 106 sc = state_color(data_matrix, markcolor)[,2];
80 } 107 }
81 rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=statecolor); 108 rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=sc);
82
83 palette(defpalette); 109 palette(defpalette);
84 dev.off(); 110 dev.off();
85 #return(statecolor);
86 } 111 }
87 112
88 # Read the input. 113 # Read the input.
89 data_matrix <- read.table(opt$input, header=T); 114 data_frame <- read.table(opt$input, comment="!", header=T);
90 create_heatmap(data_matrix); 115 create_heatmap(data_frame);