|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 build_state_color_codes_vector <- function(data_matrix, histone_mark_color, color_code_type="rgb") {
|
|
|
4 # Return vector of color code strings for each state
|
|
|
5 # in the received data_matrix. The values will be either
|
|
|
6 # rgb strings (e.g., 255,255,0) or hex code strings (e.g.,
|
|
|
7 # #FFFFFF) depending on the value of color_code_type,
|
|
|
8 # which can be one of "rgb" or "hex".
|
|
|
9 range_vector = apply(data_matrix, 1, range);
|
|
|
10 mm = NULL;
|
|
|
11 for(i in 1:dim(data_matrix)[1]) {
|
|
|
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 num_columns = dim(data_frame)[2];
|
|
|
41 num_rows = dim(data_frame)[1];
|
|
|
42 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
|
|
|
43 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
|
|
|
44 state_colors_vector = get_state_color_codes_vector(data_frame, colors=colors, color_code_type="hex");
|
|
|
45 # Open the output PDF file.
|
|
|
46 pdf(file=output_file_name);
|
|
|
47 # rownames(data_matrix) are the state indexes,
|
|
|
48 # and will look something like this:
|
|
|
49 # 0 (5.89%) 1 (91.78%) 2 (1.48%) 3 (0.86%)
|
|
|
50 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
|
|
|
51 # Set graphical parameters.
|
|
|
52 par(mar=c(6, 1, 1, 6));
|
|
|
53 # Create a vector containing the minimum and maximum values in data_matrix.
|
|
|
54 min_max_vector = range(data_matrix);
|
|
|
55 # Create a color palette.
|
|
|
56 my_palette = colorRampPalette(colors)(n=100);
|
|
|
57 default_palette = palette(my_palette);
|
|
|
58 # Plot the heatmap for the current .para / .state combination.
|
|
|
59 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);
|
|
|
60 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
|
|
|
61 axis(4, at=1:num_rows-0.5, labels=rownames(data_matrix), las=2);
|
|
|
62 col = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100);
|
|
|
63 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);
|
|
|
64 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);
|
|
|
65 palette(default_palette);
|
|
|
66 dev.off();
|
|
|
67 }
|
|
|
68
|
|
|
69 get_state_color_codes_vector <- function(data_frame, colors=c("white", "dark blue"), color_code_type="rgb") {
|
|
|
70 # Return a vector of color strings for each row in data_frame.
|
|
|
71 # These string will either be rgb (e.g., 255,255,0) or hex codes
|
|
|
72 # (e.g., #FFFFFF), depending on the value of color_code_type.
|
|
|
73 num_columns = dim(data_frame)[2];
|
|
|
74 num_rows = dim(data_frame)[1];
|
|
|
75 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
|
|
|
76 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
|
|
|
77 # colnames(data_matrix) will look something like this:
|
|
|
78 # H3K4me3 H3K4me1 DNase H3K79me2
|
|
|
79 colnames(data_matrix) = colnames(data_frame)[1+1:p];
|
|
|
80 histone_marks = colnames(data_matrix);
|
|
|
81 histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
|
|
|
82 # Specify colors for common feature names like "h3k4me3".
|
|
|
83 # These are histone marks frequently used to identify
|
|
|
84 # promoter activities in a cell, and are often displayed
|
|
|
85 # in shades of red.
|
|
|
86 for(i in 1:length(histone_marks)) {
|
|
|
87 if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
|
|
|
88 histone_mark_color[i,] = c(255, 0, 0);
|
|
|
89 }
|
|
|
90 if(regexpr("h3k4me2", tolower(histone_marks[i])) > 0) {
|
|
|
91 histone_mark_color[i,] = c(250, 100, 0);
|
|
|
92 }
|
|
|
93 if(regexpr("h3k4me1", tolower(histone_marks[i])) > 0) {
|
|
|
94 histone_mark_color[i,] = c(250, 250, 0);
|
|
|
95 }
|
|
|
96 if(regexpr("h3k36me3", tolower(histone_marks[i]))>0) {
|
|
|
97 histone_mark_color[i,] = c(0, 150, 0);
|
|
|
98 }
|
|
|
99 if(regexpr("h2a", tolower(histone_marks[i])) > 0) {
|
|
|
100 histone_mark_color[i,] = c(0, 150, 150);
|
|
|
101 }
|
|
|
102 if(regexpr("dnase", tolower(histone_marks[i])) > 0) {
|
|
|
103 histone_mark_color[i,] = c(0, 200, 200);
|
|
|
104 }
|
|
|
105 if(regexpr("h3k9ac", tolower(histone_marks[i])) > 0) {
|
|
|
106 histone_mark_color[i,] = c(250, 0, 200);
|
|
|
107 }
|
|
|
108 if(regexpr("h3k9me3", tolower(histone_marks[i])) > 0) {
|
|
|
109 histone_mark_color[i,] = c(100, 100, 100);
|
|
|
110 }
|
|
|
111 if(regexpr("h3k27ac", tolower(histone_marks[i])) > 0) {
|
|
|
112 histone_mark_color[i,] = c(250, 150, 0);
|
|
|
113 }
|
|
|
114 if(regexpr("h3k27me3", tolower(histone_marks[i])) > 0) {
|
|
|
115 histone_mark_color[i,] = c(0, 0, 200);
|
|
|
116 }
|
|
|
117 if(regexpr("h3k79me2", tolower(histone_marks[i])) > 0) {
|
|
|
118 histone_mark_color[i,] = c(200, 0, 200);
|
|
|
119 }
|
|
|
120 if(regexpr("h4k20me1", tolower(histone_marks[i])) > 0) {
|
|
|
121 histone_mark_color[i,] = c(50, 200, 50);
|
|
|
122 }
|
|
|
123 if(regexpr("ctcf", tolower(histone_marks[i])) > 0) {
|
|
|
124 histone_mark_color[i,] = c(200, 0, 250);
|
|
|
125 }
|
|
|
126 state_colors_vector = build_state_color_codes_vector(data_matrix, histone_mark_color, color_code_type=color_code_type);
|
|
|
127 }
|
|
|
128 return(state_colors_vector);
|
|
|
129 } |
