Mercurial > repos > greg > ideas_genome_tracks

comparison create_heatmap.R @ 66:f0c488b73dbc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Wed, 20 Dec 2017 11:26:04 -0500
parents
children efbc9e58c4eb
comparison
equal deleted inserted replaced
65:794d2104f83d 66:f0c488b73dbc
1 #!/usr/bin/env Rscript
2
3 create_heatmap<-function(data_frame, output_file_name=NULL) {
4 # Plot a heatmap for a .para / .state combination
5 # based on the received data_frame which was created
6 # by reading the .para file.
7 num_columns = dim(data_frame)[2];
8 num_rows = dim(data_frame)[1];
9 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
10 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
11 colnames(data_matrix) = colnames(data_frame)[1+1:p];
12 histone_marks = colnames(data_matrix);
13 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.
19 par(mar=c(6, 1, 1, 6));
20 # Create a vector containing the minimum and maximum values in data_matrix.
21 min_max_vector = range(data_matrix);
22 # Create a color palette.
23 my_palette = colorRampPalette(c("white", "dark blue"))(n=100);
24 defpalette = palette(my_palette);
25 # 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);
27 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);
29 color = 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=color);
31 histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
32
33 # Specify a color for common feature names like "h3k4me3".
34 # These are histone marks frequently used to identify
35 # promoter activities in a cell, and are often displayed
36 # in shades of red.
37 for(i in 1:length(histone_marks)) {
38 if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
39 histone_mark_color[i,] = c(255, 0, 0);
40 }
41 if(regexpr("h3k4me2", tolower(histone_marks[i])) > 0) {
42 histone_mark_color[i,] = c(250, 100, 0);
43 }
44 if(regexpr("h3k4me1", tolower(histone_marks[i])) > 0) {
45 histone_mark_color[i,] = c(250, 250, 0);
46 }
47 if(regexpr("h3k36me3", tolower(histone_marks[i]))>0) {
48 histone_mark_color[i,] = c(0, 150, 0);
49 }
50 if(regexpr("h2a", tolower(histone_marks[i])) > 0) {
51 histone_mark_color[i,] = c(0, 150, 150);
52 }
53 if(regexpr("dnase", tolower(histone_marks[i])) > 0) {
54 histone_mark_color[i,] = c(0, 200, 200);
55 }
56 if(regexpr("h3k9ac", tolower(histone_marks[i])) > 0) {
57 histone_mark_color[i,] = c(250, 0, 200);
58 }
59 if(regexpr("h3k9me3", tolower(histone_marks[i])) > 0) {
60 histone_mark_color[i,] = c(100, 100, 100);
61 }
62 if(regexpr("h3k27ac", tolower(histone_marks[i])) > 0) {
63 histone_mark_color[i,] = c(250, 150, 0);
64 }
65 if(regexpr("h3k27me3", tolower(histone_marks[i])) > 0) {
66 histone_mark_color[i,] = c(0, 0, 200);
67 }
68 if(regexpr("h3k79me2", tolower(histone_marks[i])) > 0) {
69 histone_mark_color[i,] = c(200, 0, 200);
70 }
71 if(regexpr("h4k20me1", tolower(histone_marks[i])) > 0) {
72 histone_mark_color[i,] = c(50, 200, 50);
73 }
74 if(regexpr("ctcf", tolower(histone_marks[i])) > 0) {
75 histone_mark_color[i,] = c(200, 0, 250);
76 }
77 state_color = get_state_color(data_matrix, histone_mark_color)[,2];
78 }
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);
80 palette(defpalette);
81 if (!is.null(output_file_name)) {
82 dev.off();
83 }
84 return(state_color);
85 }
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 }