Mercurial > repos > greg > ideas

comparison create_heatmap.R @ 133:7d49fc24e06a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Fri, 15 Dec 2017 10:25:34 -0500
parents 7229e52fa8e1
children fc94a1ce21eb
comparison
equal deleted inserted replaced
132:8ce93420010c 133:7d49fc24e06a
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 state_color <- function(statemean, markcolor=NULL) 14 create_heatmap<-function(data_frame, output_file_name) {
15 { 15 # Plot a heatmap for a .para / .state combination
16 if(length(markcolor) == 0) { 16 # based on the received data_frame which was created
17 markcolor = rep("", dim(statemean)[2]); 17 # by reading the .para file.
18 markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:dim(statemean)[2]-1) / dim(statemean)[2], 1, 1); 18 num_columns = dim(data_frame)[2];
19 markcolor = t(col2rgb(markcolor)); 19 num_rows = dim(data_frame)[1];
20 p = (sqrt(9 + 8 * (num_columns - 1)) - 3) / 2;
21 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
22 column_names(data_matrix) = column_names(data_frame)[1+1:p];
23 marks = column_names(data_matrix);
24 row_names(data_matrix) = paste(1:num_rows, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
25
26 # Open the output PDF file.
27 pdf(file=output_file_name);
28 # Set graphical parameters.
29 par(mar=c(6, 1, 1, 6));
30 # Create a vector containing the minimum and maximum values in data_matrix.
31 min_max_vector = range(data_matrix);
32 # Define colors for the palette.
33 colors = 0:100 / 100 * (min_max_vector[2] - min_max_vector[1]) + min_max_vector[1];
34 # Create the color palette.
35 my_palette = colorRampPalette(c("white", "dark blue"))(n=100);
36 defpalette = palette(my_palette);
37 # Plot the heatmap for the current .para / .state combination.
38 plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
39 axis(1, at=1:p-0.5, labels=column_names(data_matrix), las=2);
40 axis(4, at=1:num_rows-0.5, labels=row_names(data_matrix), las=2);
41 rect(rep(1:p-1, num_rows), rep(1:num_rows, each=p), rep(1:p, num_rows), rep(1:num_rows, each=p),
42 col=round((t(data_matrix)-min_max_vector[1])/(min_max_vector[2]-min_max_vector[1])*100));
43 markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
44
45 for(i in 1:length(marks)) {
46 if(regexpr("h3k4me3", tolower(marks[i])) > 0) {
47 markcolor[i,] = c(255, 0, 0);
48 }
49 if(regexpr("h3k4me2", tolower(marks[i])) > 0) {
50 markcolor[i,] = c(250, 100, 0);
51 }
52 if(regexpr("h3k4me1", tolower(marks[i])) > 0) {
53 markcolor[i,] = c(250, 250, 0);
54 }
55 if(regexpr("h3k36me3", tolower(marks[i]))>0) {
56 markcolor[i,] = c(0, 150, 0);
57 }
58 if(regexpr("h2a", tolower(marks[i])) > 0) {
59 markcolor[i,] = c(0, 150, 150);
60 }
61 if(regexpr("dnase", tolower(marks[i])) > 0) {
62 markcolor[i,] = c(0, 200, 200);
63 }
64 if(regexpr("h3k9ac", tolower(marks[i])) > 0) {
65 markcolor[i,] = c(250, 0, 200);
66 }
67 if(regexpr("h3k9me3", tolower(marks[i])) > 0) {
68 markcolor[i,] = c(100, 100, 100);
69 }
70 if(regexpr("h3k27ac", tolower(marks[i])) > 0) {
71 markcolor[i,] = c(250, 150, 0);
72 }
73 if(regexpr("h3k27me3", tolower(marks[i])) > 0) {
74 markcolor[i,] = c(0, 0, 200);
75 }
76 if(regexpr("h3k79me2", tolower(marks[i])) > 0) {
77 markcolor[i,] = c(200, 0, 200);
78 }
79 if(regexpr("h4k20me1", tolower(marks[i])) > 0) {
80 markcolor[i,] = c(50, 200, 50);
81 }
82 if(regexpr("ctcf", tolower(marks[i])) > 0) {
83 markcolor[i,] = c(200, 0, 250);
84 }
85 state_color = get_state_color(data_matrix, markcolor)[,2];
20 } 86 }
87 rect(rep(p+0.2,num_rows), 1:num_rows-0.8, rep(p+0.8,num_rows), 1:num_rows-0.2, col=state_color);
88 palette(defpalette);
89 dev.off();
90 }
21 91
22 rg = apply(statemean, 1, range); 92 get_state_color<-function(statemean, markcolor) {
23 mm = NULL; 93 rg=apply(statemean, 1, range);
94 mm=NULL;
24 for(i in 1:dim(statemean)[1]) { 95 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)); 96 mm = rbind(mm, (statemean[i,]-rg[1,i]+1e-10)/(rg[2,i]-rg[1,i]+1e-10));
26 } 97 }
27 mm = mm^6; 98 mm = mm^5;
28 if(dim(mm)[2] > 1) { 99 if(dim(mm)[2] > 1) {
29 mm = mm / (apply(mm, 1, sum) + 1e-10); 100 mm = mm / (apply(mm, 1, sum) + 1e-10);
30 } 101 }
31 mycol = mm%*%markcolor; 102 mycol = mm%*%markcolor;
32 s = apply(statemean, 1, max); 103 s = apply(statemean, 1, max);
33 s = (s - min(s)) / (max(s) - min(s) + 1e-10); 104 s = (s-min(s)) / (max(s)-min(s) + 1e-10);
105 mycol = round(255-(255-mycol)*s/0.5);
106 mycol[mycol<0] = 0;
107 rt = paste(mycol[,1], mycol[,2], mycol[,3], sep=",");
34 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])})); 108 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])}); 109 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); 110 rt = cbind(rt, h);
38 return(rt); 111 return(rt);
39 }
40
41 create_heatmap<-function(data_frame, output_file_name, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) {
42 k = dim(data_frame)[2];
43 l = dim(data_frame)[1];
44 p = (sqrt(9 + 8 * (k - 1)) - 3) / 2;
45 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
46 colnames(data_matrix) = colnames(data_frame)[1+1:p];
47 marks = colnames(data_matrix);
48 rownames(data_matrix) = paste(1:l," (", round(data_frame[,1] / sum(data_frame[,1]) * 10000) / 100, "%)", sep="");
49 pdf(file=output_file_name);
50 par(mar=c(6, 1, 1, 6));
51 rg = range(data_matrix);
52 colors = 0:100 / 100 * (rg[2] - rg[1]) + rg[1];
53 my_palette = colorRampPalette(cols)(n=100);
54 defpalette = palette(my_palette);
55
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);
57 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
58 axis(4, at=1:l-0.5, labels=rownames(data_matrix), las=2);
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));
60
61 if(length(statecolor)==0) {
62 if(length(markcolor)==0) {
63 markcolor=t(col2rgb(terrain.colors(ceiling(p))[1:p]));
64 for(i in 1:length(marks)) {
65 if(regexpr("h3k4me3",tolower(marks[i]))>0) {
66 markcolor[i,]=c(255,0,0);
67 }
68 if(regexpr("h3k4me2",tolower(marks[i]))>0) {
69 markcolor[i,]=c(250,100,0);
70 }
71 if(regexpr("h3k4me1",tolower(marks[i]))>0) {
72 markcolor[i,]=c(250,250,0);
73 }
74 if(regexpr("h3k36me3",tolower(marks[i]))>0) {
75 markcolor[i,]=c(0,150,0);
76 }
77 if(regexpr("h2a",tolower(marks[i]))>0) {
78 markcolor[i,]=c(0,150,150);
79 }
80 if(regexpr("dnase",tolower(marks[i]))>0) {
81 markcolor[i,]=c(0,200,200);
82 }
83 if(regexpr("h3k9ac",tolower(marks[i]))>0) {
84 markcolor[i,]=c(250,0,200);
85 }
86 if(regexpr("h3k9me3",tolower(marks[i]))>0) {
87 markcolor[i,]=c(100,100,100);
88 }
89 if(regexpr("h3k27ac",tolower(marks[i]))>0) {
90 markcolor[i,]=c(250,150,0);
91 }
92 if(regexpr("h3k27me3",tolower(marks[i]))>0) {
93 markcolor[i,]=c(0,0,200);
94 }
95 if(regexpr("h3k79me2",tolower(marks[i]))>0) {
96 markcolor[i,]=c(200,0,200);
97 }
98 if(regexpr("h4k20me1",tolower(marks[i]))>0) {
99 markcolor[i,]=c(50,200,50);
100 }
101 if(regexpr("ctcf",tolower(marks[i]))>0) {
102 markcolor[i,]=c(200,0,250);
103 }
104 }
105 }
106 sc = state_color(data_matrix, markcolor)[,2];
107 }
108 rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=sc);
109 palette(defpalette);
110 dev.off();
111 } 112 }
112 113
113 # Read the inputs. 114 # Read the inputs.
114 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE); 115 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE);
115 for (i in 1:length(para_files)) { 116 for (i in 1:length(para_files)) {
116 para_file <- para_files[i]; 117 para_file <- para_files[i];
117 para_file_base_name <- strsplit(para_file, split="/")[[1]][2] 118 para_file_base_name <- strsplit(para_file, split="/")[[1]][2]
118 output_file_name <- gsub(".para", ".pdf", para_file_base_name) 119 output_file_name <- gsub(".para", ".pdf", para_file_base_name)
119 output_file_path <- paste(opt$output_dir, "/", output_file_name, sep=""); 120 output_file_path <- paste(opt$output_dir, output_file_name, sep="/");
120 data_frame <- read.table(para_file, comment="!", header=T); 121 data_frame <- read.table(para_file, comment="!", header=T);
121 create_heatmap(data_frame, output_file_path); 122 create_heatmap(data_frame, output_file_path);
122 } 123 }