Mercurial > repos > greg > ideas

comparison create_heatmap.R @ 147:387b460ddd43 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 11 Jan 2018 10:21:13 -0500
parents 9e0b4ceba74a
children 0ba72d5ca209
comparison
equal deleted inserted replaced
146:9e0b4ceba74a 147:387b460ddd43
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 create_heatmap<-function(data_frame, output_file_name=NULL) { 3 build_state_color_codes_vector <- function(data_matrix, histone_mark_color, color_code_type="rgb") {
4 # Plot a heatmap for a .para / .state combination 4 # Return vector of color code strings for each state
5 # based on the received data_frame which was created 5 # in the received data_matrix. The values will be either
6 # by reading the .para file. 6 # rgb strings (e.g., 255,255,0) or hex code strings (e.g.,
7 num_columns = dim(data_frame)[2]; 7 # #FFFFFF) depending on the value of color_code_type,
8 num_rows = dim(data_frame)[1]; 8 # which can be one of "rgb" or "hex".
9 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2; 9 range_vector = apply(data_matrix, 1, range);
10 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]); 10 mm = NULL;
11 colnames(data_matrix) = colnames(data_frame)[1+1:p]; 11 for(i in 1:dim(data_matrix)[1]) {
12 histone_marks = colnames(data_matrix); 12 range_val1 = range_vector[1, i] + 1e-10;
13 range_val2 = range_vector[2, i];
14 mm = rbind(mm, (data_matrix[i,] - range_val1) / (range_val2 - range_val1));
15 }
16 mm = mm^5;
17 if(dim(mm)[2] > 1) {
18 mm = mm / (apply(mm, 1, sum) + 1e-10);
19 }
20 state_color = mm%*%histone_mark_color;
21 s = apply(data_matrix, 1, max);
22 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
23 state_color = round(255 - (255 - state_color) * s/0.5);
24 state_color[state_color<0] = 0;
25 if (identical(color_code_type, "rgb")) {
26 # Here rgb_values is something like 255,255,255 217,98,0.
27 state_colors_vector = paste(state_color[,1], state_color[,2], state_color[,3], sep=",");
28 } else {
29 # Here hex_code_strings is something like #FFFFFF #D96200
30 # which is a one-to-one map to the above rgb_values.
31 hex_code_strings = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
32 state_colors_vector = apply(hex_code_strings, 1, function(x){hsv(x[1], x[2], x[3])});
33 }
34 return(state_colors_vector);
35 }
36
37 create_heatmap <- function(data_frame, output_file_name, colors=c("white", "dark blue")) {
38 # Plot a heatmap for a .para / .state combination based on the
39 # received data_frame which was created by reading the .para file.
40 state_colors_vector = get_state_color_codes_vector(data_frame, colors=colors, color_code_type="hex");
41 # Open the output PDF file.
42 pdf(file=output_file_name);
43 # rownames(data_matrix) are the state indexes,
44 # and will look something like this:
45 # 0 (5.89%) 1 (91.78%) 2 (1.48%) 3 (0.86%)
13 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep=""); 46 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
14 if (!is.null(output_file_name)) {
15 # Open the output PDF file.
16 pdf(file=output_file_name);
17 }
18 # Set graphical parameters. 47 # Set graphical parameters.
19 par(mar=c(6, 1, 1, 6)); 48 par(mar=c(6, 1, 1, 6));
20 # Create a vector containing the minimum and maximum values in data_matrix. 49 # Create a vector containing the minimum and maximum values in data_matrix.
21 min_max_vector = range(data_matrix); 50 min_max_vector = range(data_matrix);
22 # Create a color palette. 51 # Create a color palette.
23 my_palette = colorRampPalette(c("white", "dark blue"))(n=100); 52 my_palette = colorRampPalette(colors)(n=100);
24 defpalette = palette(my_palette); 53 default_palette = palette(my_palette);
25 # Plot the heatmap for the current .para / .state combination. 54 # Plot the heatmap for the current .para / .state combination.
26 plot(NA, NA, xlim=c(0, p+0.7), ylim=c(0, num_rows), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F); 55 plot(NA, NA, xlim=c(0, p+0.7), ylim=c(0, num_rows), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
27 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2); 56 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
28 axis(4, at=1:num_rows-0.5, labels=rownames(data_matrix), las=2); 57 axis(4, at=1:num_rows-0.5, labels=rownames(data_matrix), las=2);
29 col = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100); 58 col = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100);
30 rect(rep(1:p-1, num_rows), rep(1:num_rows-1, each=p), rep(1:p, num_rows), rep(1:num_rows, each=p), col=col); 59 rect(rep(1:p-1, num_rows), rep(1:num_rows-1, each=p), rep(1:p, num_rows), rep(1:num_rows, each=p), col=col);
60 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_colors_vector);
61 palette(default_palette);
62 dev.off();
63 }
64
65 get_state_color_codes_vector <- function(data_frame, colors=c("white", "dark blue"), color_code_type="rgb") {
66 # Return a vector of color strings for each row in data_frame.
67 # These string will either be rgb (e.g., 255,255,0) or hex codes
68 # (e.g., #FFFFFF), depending on the value of color_code_type.
69 num_columns = dim(data_frame)[2];
70 num_rows = dim(data_frame)[1];
71 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
72 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
73 # colnames(data_matrix) will look something like this:
74 # H3K4me3 H3K4me1 DNase H3K79me2
75 colnames(data_matrix) = colnames(data_frame)[1+1:p];
76 histone_marks = colnames(data_matrix);
31 histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p])); 77 histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
32 78 # Specify colors for common feature names like "h3k4me3".
33 # Specify a color for common feature names like "h3k4me3".
34 # These are histone marks frequently used to identify 79 # These are histone marks frequently used to identify
35 # promoter activities in a cell, and are often displayed 80 # promoter activities in a cell, and are often displayed
36 # in shades of red. 81 # in shades of red.
37 for(i in 1:length(histone_marks)) { 82 for(i in 1:length(histone_marks)) {
38 if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) { 83 if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
72 histone_mark_color[i,] = c(50, 200, 50); 117 histone_mark_color[i,] = c(50, 200, 50);
73 } 118 }
74 if(regexpr("ctcf", tolower(histone_marks[i])) > 0) { 119 if(regexpr("ctcf", tolower(histone_marks[i])) > 0) {
75 histone_mark_color[i,] = c(200, 0, 250); 120 histone_mark_color[i,] = c(200, 0, 250);
76 } 121 }
77 state_color = get_state_color(data_matrix, histone_mark_color)[,]; 122 state_colors_vector = build_state_color_codes_vector(data_matrix, histone_mark_color, color_code_type=color_code_type);
78 } 123 }
79 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[,2]); 124 return(state_colors_vector);
80 palette(defpalette);
81 if (!is.null(output_file_name)) {
82 dev.off();
83 }
84 return(state_color);
85 } 125 }
86
87 get_state_color <- function(data_matrix, histone_mark_color) {
88 range_vector = apply(data_matrix, 1, range);
89 mm = NULL;
90 for(i in 1:dim(data_matrix)[1]) {
91 range_val1 = range_vector[1, i] + 1e-10;
92 range_val2 = range_vector[2, i];
93 mm = rbind(mm, (data_matrix[i,] - range_val1) / (range_val2 - range_val1));
94 }
95 mm = mm^5;
96 if(dim(mm)[2] > 1) {
97 mm = mm / (apply(mm, 1, sum) + 1e-10);
98 }
99 state_color = mm%*%histone_mark_color;
100 s = apply(data_matrix, 1, max);
101 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
102 state_color = round(255 - (255 - state_color) * s/0.5);
103 state_color[state_color<0] = 0;
104 rt = paste(state_color[,1], state_color[,2], state_color[,3], sep=",");
105 h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
106 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
107 rt = cbind(rt, h);
108 return(rt);
109 }