Mercurial > repos > greg > ideas

comparison create_heatmap.R @ 136:8d41990f4fc5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 18 Dec 2017 14:30:14 -0500
parents fc94a1ce21eb
children 3180393013cc
comparison
equal deleted inserted replaced
135:fc94a1ce21eb 136:8d41990f4fc5
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-i", "--input_dir"), action="store", dest="input_dir", help="IDEAS para files directory"), 6 make_option(c("-i", "--input_dir"), action="store", dest="input_dir", help="IDEAS para files directory"),
7 make_option(c("-o", "--output_dir"), action="store", dest="output_dir", help="PDF output directory") 7 make_option(c("-o", "--output_dir"), action="store", dest="output_dir", help="PDF output directory")
8 ) 8 )
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
14 create_heatmap<-function(data_frame, output_file_name) { 14 create_heatmap<-function(data_frame, output_file_name) {
15 # Plot a heatmap for a .para / .state combination 15 # Plot a heatmap for a .para / .state combination
16 # based on the received data_frame which was created 16 # based on the received data_frame which was created
17 # by reading the .para file. 17 # by reading the .para file.
18 num_columns = dim(data_frame)[2]; 18 num_columns = dim(data_frame)[2];
19 cat("num_columns: ", num_columns, "\n");
19 num_rows = dim(data_frame)[1]; 20 num_rows = dim(data_frame)[1];
20 p = (sqrt(9 + 8 * (num_columns - 1)) - 3) / 2; 21 cat("num_rows: ", num_rows, "\n");
22 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
23 cat("9 + 8 * (num_columns-1): ", 9 + 8 * (num_columns-1), "\n");
24 cat("sqrt(9 + 8 * (num_columns-1)): ", sqrt(9 + 8 * (num_columns-1)), "\n");
25 cat("p: ", p, "\n");
21 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]); 26 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
27 cat("dim(data_matrix)[1]: ", dim(data_matrix)[1], "\n");
28 cat("dim(data_matrix)[2]: ", dim(data_matrix)[2], "\n");
22 colnames(data_matrix) = colnames(data_frame)[1+1:p]; 29 colnames(data_matrix) = colnames(data_frame)[1+1:p];
23 marks = colnames(data_matrix); 30 cat("colnames(data_matrix): ", colnames(data_matrix), "\n");
24 rownames(data_matrix) = paste(1:num_rows, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep=""); 31 histone_marks = colnames(data_matrix);
25 32 cat("histone_marks: ", histone_marks, "\n");
33 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
34 cat("rownames(data_matrix): ", rownames(data_matrix), "\n");
26 # Open the output PDF file. 35 # Open the output PDF file.
27 pdf(file=output_file_name); 36 pdf(file=output_file_name);
28 # Set graphical parameters. 37 # Set graphical parameters.
29 par(mar=c(6, 1, 1, 6)); 38 par(mar=c(6, 1, 1, 6));
30 # Create a vector containing the minimum and maximum values in data_matrix. 39 # Create a vector containing the minimum and maximum values in data_matrix.
31 min_max_vector = range(data_matrix); 40 min_max_vector = range(data_matrix);
32 # Define colors for the palette. 41 cat("min_max_vector: ", min_max_vector, "\n");
33 colors = 0:100 / 100 * (min_max_vector[2] - min_max_vector[1]) + min_max_vector[1]; 42 # Create a color palette.
34 # Create the color palette.
35 my_palette = colorRampPalette(c("white", "dark blue"))(n=100); 43 my_palette = colorRampPalette(c("white", "dark blue"))(n=100);
36 defpalette = palette(my_palette); 44 defpalette = palette(my_palette);
37 # Plot the heatmap for the current .para / .state combination. 45 # Plot the heatmap for the current .para / .state combination.
38 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); 46 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);
39 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2); 47 axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
40 axis(4, at=1:num_rows-0.5, labels=rownames(data_matrix), las=2); 48 axis(4, at=1:num_rows-0.5, labels=rownames(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), 49 color = round((t(data_matrix) - min_max_vector[1]) / (min_max_vector[2] - min_max_vector[1]) * 100);
42 col=round((t(data_matrix)-min_max_vector[1])/(min_max_vector[2]-min_max_vector[1])*100)); 50 cat("color: ", color, "\n");
43 markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p])); 51 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);
52 histone_mark_color = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
53 cat("histone_mark_color: ", histone_mark_color, "\n");
44 54
45 for(i in 1:length(marks)) { 55 # Specify a color for common feature names like "h3k4me3".
46 if(regexpr("h3k4me3", tolower(marks[i])) > 0) { 56 # These are histone marks frequently used to identify
47 markcolor[i,] = c(255, 0, 0); 57 # promoter activities in a cell, and are often displayed
58 # in shades of red.
59 for(i in 1:length(histone_marks)) {
60 cat("histone_marks[i]: ", histone_marks[i], "\n");
61 if(regexpr("h3k4me3", tolower(histone_marks[i])) > 0) {
62 histone_mark_color[i,] = c(255, 0, 0);
63 cat("histone_mark_color[i]: ", histone_mark_color[i], "\n");
48 } 64 }
49 if(regexpr("h3k4me2", tolower(marks[i])) > 0) { 65 if(regexpr("h3k4me2", tolower(histone_marks[i])) > 0) {
50 markcolor[i,] = c(250, 100, 0); 66 histone_mark_color[i,] = c(250, 100, 0);
51 } 67 }
52 if(regexpr("h3k4me1", tolower(marks[i])) > 0) { 68 if(regexpr("h3k4me1", tolower(histone_marks[i])) > 0) {
53 markcolor[i,] = c(250, 250, 0); 69 histone_mark_color[i,] = c(250, 250, 0);
54 } 70 }
55 if(regexpr("h3k36me3", tolower(marks[i]))>0) { 71 if(regexpr("h3k36me3", tolower(histone_marks[i]))>0) {
56 markcolor[i,] = c(0, 150, 0); 72 histone_mark_color[i,] = c(0, 150, 0);
57 } 73 }
58 if(regexpr("h2a", tolower(marks[i])) > 0) { 74 if(regexpr("h2a", tolower(histone_marks[i])) > 0) {
59 markcolor[i,] = c(0, 150, 150); 75 histone_mark_color[i,] = c(0, 150, 150);
60 } 76 }
61 if(regexpr("dnase", tolower(marks[i])) > 0) { 77 if(regexpr("dnase", tolower(histone_marks[i])) > 0) {
62 markcolor[i,] = c(0, 200, 200); 78 histone_mark_color[i,] = c(0, 200, 200);
63 } 79 }
64 if(regexpr("h3k9ac", tolower(marks[i])) > 0) { 80 if(regexpr("h3k9ac", tolower(histone_marks[i])) > 0) {
65 markcolor[i,] = c(250, 0, 200); 81 histone_mark_color[i,] = c(250, 0, 200);
66 } 82 }
67 if(regexpr("h3k9me3", tolower(marks[i])) > 0) { 83 if(regexpr("h3k9me3", tolower(histone_marks[i])) > 0) {
68 markcolor[i,] = c(100, 100, 100); 84 histone_mark_color[i,] = c(100, 100, 100);
69 } 85 }
70 if(regexpr("h3k27ac", tolower(marks[i])) > 0) { 86 if(regexpr("h3k27ac", tolower(histone_marks[i])) > 0) {
71 markcolor[i,] = c(250, 150, 0); 87 histone_mark_color[i,] = c(250, 150, 0);
72 } 88 }
73 if(regexpr("h3k27me3", tolower(marks[i])) > 0) { 89 if(regexpr("h3k27me3", tolower(histone_marks[i])) > 0) {
74 markcolor[i,] = c(0, 0, 200); 90 histone_mark_color[i,] = c(0, 0, 200);
75 } 91 }
76 if(regexpr("h3k79me2", tolower(marks[i])) > 0) { 92 if(regexpr("h3k79me2", tolower(histone_marks[i])) > 0) {
77 markcolor[i,] = c(200, 0, 200); 93 histone_mark_color[i,] = c(200, 0, 200);
78 } 94 }
79 if(regexpr("h4k20me1", tolower(marks[i])) > 0) { 95 if(regexpr("h4k20me1", tolower(histone_marks[i])) > 0) {
80 markcolor[i,] = c(50, 200, 50); 96 histone_mark_color[i,] = c(50, 200, 50);
81 } 97 }
82 if(regexpr("ctcf", tolower(marks[i])) > 0) { 98 if(regexpr("ctcf", tolower(histone_marks[i])) > 0) {
83 markcolor[i,] = c(200, 0, 250); 99 histone_mark_color[i,] = c(200, 0, 250);
84 } 100 }
85 state_color = get_state_color(data_matrix, markcolor)[,2]; 101 state_color = get_state_color(data_matrix, histone_mark_color)[,2];
102 cat("state_color: ", state_color, "\n");
86 } 103 }
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); 104 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); 105 palette(defpalette);
89 dev.off(); 106 dev.off();
90 } 107 }
91 108
92 get_state_color<-function(statemean, markcolor) { 109 get_state_color <- function(data_matrix, histone_mark_color) {
93 rg=apply(statemean, 1, range); 110 cat("XXX histone_mark_color: ", histone_mark_color, "\n");
94 mm=NULL; 111 #histone_mark_color = rep("", dim(data_matrix)[2]);
95 for(i in 1:dim(statemean)[1]) { 112 #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
96 mm = rbind(mm, (statemean[i,]-rg[1,i]+1e-10)/(rg[2,i]-rg[1,i]+1e-10)); 113 #cat("XXX order(apply(data_matrix, 2, sd), decreasing=T): ", order(apply(data_matrix, 2, sd), decreasing=T), "\n");
114 #cat("XXX dim(data_matrix)[2]-1: ", dim(data_matrix)[2]-1, "\n");
115 #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
116 #histone_mark_color[order(apply(data_matrix, 2, sd), decreasing=T)] = hsv((1:dim(data_matrix)[2]-1) / dim(data_matrix)[2], 1, 1);
117 #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
118 #histone_mark_color = t(col2rgb(histone_mark_color));
119 #cat("XXX histone_mark_color: ", histone_mark_color, "\n");
120 rg = apply(data_matrix, 1, range);
121 cat("XXX rg: ", rg, "\n");
122 mm = NULL;
123 cat("XXX dim(data_matrix): ", dim(data_matrix), "\n");
124 cat("XXX dim(data_matrix)[1]: ", dim(data_matrix)[1], "\n");
125 for(i in 1:dim(data_matrix)[1]) {
126 mm = rbind(mm, (data_matrix[i,] - rg[1, i] + 1e-10) / (rg[2, i] -rg[1, i] + 1e-10));
127 cat("XXX mm: ", mm, "\n");
97 } 128 }
98 mm = mm^5; 129 mm = mm^5;
130 cat("XXX dim(mm): ", dim(mm), "\n");
131 cat("XXX dim(mm)[2]: ", dim(mm)[2], "\n");
99 if(dim(mm)[2] > 1) { 132 if(dim(mm)[2] > 1) {
100 mm = mm / (apply(mm, 1, sum) + 1e-10); 133 mm = mm / (apply(mm, 1, sum) + 1e-10);
134 cat("XXX mm: ", mm, "\n");
101 } 135 }
102 mycol = mm%*%markcolor; 136 state_color = mm%*%histone_mark_color;
103 s = apply(statemean, 1, max); 137 cat("XXX state_color: ", state_color, "\n");
104 s = (s-min(s)) / (max(s)-min(s) + 1e-10); 138 s = apply(data_matrix, 1, max);
105 mycol = round(255-(255-mycol)*s/0.5); 139 cat("XXX s: ", s, "\n");
106 mycol[mycol<0] = 0; 140 s = (s - min(s)) / (max(s) - min(s) + 1e-10);
107 rt = paste(mycol[,1], mycol[,2], mycol[,3], sep=","); 141 cat("XXX s: ", s, "\n");
108 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])})); 142 state_color = round(255 - (255 - state_color) * s/0.5);
143 cat("XXX state_color: ", state_color, "\n");
144 state_color[state_color<0] = 0;
145 cat("XXX state_color: ", state_color, "\n");
146 rt = paste(state_color[,1], state_color[,2], state_color[,3], sep=",");
147 cat("XXX rt: ", rt, "\n");
148 h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
149 cat("XXX h: ", h, "\n");
109 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])}); 150 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
151 cat("XXX h: ", h, "\n");
110 rt = cbind(rt, h); 152 rt = cbind(rt, h);
153 cat("XXX rt: ", rt, "\n");
111 return(rt); 154 return(rt);
112 } 155 }
113 156
114 # Read the inputs. 157 # Read the inputs.
115 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE); 158 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE);