Mercurial > repos > greg > ideas

comparison create_heatmap.R @ 141:a976dd6fcd1b 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:24:23 -0500
parents 17e8829bbae2
children 8cd08d0aeb69
comparison
equal deleted inserted replaced
140:17e8829bbae2 141:a976dd6fcd1b
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 create_heatmap<-function(data_frame, output_file_name=NULL) {
4
5 option_list <- list(
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")
8 )
9
10 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
11 args <- parse_args(parser, positional_arguments=TRUE)
12 opt <- args$options
13
14 create_heatmap<-function(data_frame, output_file_name) {
15 # Plot a heatmap for a .para / .state combination 4 # Plot a heatmap for a .para / .state combination
16 # based on the received data_frame which was created 5 # based on the received data_frame which was created
17 # by reading the .para file. 6 # by reading the .para file.
18 num_columns = dim(data_frame)[2]; 7 num_columns = dim(data_frame)[2];
19 num_rows = dim(data_frame)[1]; 8 num_rows = dim(data_frame)[1];
20 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2; 9 p = (sqrt(9 + 8 * (num_columns-1)) - 3) / 2;
21 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]); 10 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
22 colnames(data_matrix) = colnames(data_frame)[1+1:p]; 11 colnames(data_matrix) = colnames(data_frame)[1+1:p];
23 histone_marks = colnames(data_matrix); 12 histone_marks = colnames(data_matrix);
24 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep=""); 13 rownames(data_matrix) = paste(1:num_rows-1, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
25 # Open the output PDF file. 14 if (!is.null(output_file_name)) {
26 pdf(file=output_file_name); 15 # Open the output PDF file.
16 pdf(file=output_file_name);
17 }
27 # Set graphical parameters. 18 # Set graphical parameters.
28 par(mar=c(6, 1, 1, 6)); 19 par(mar=c(6, 1, 1, 6));
29 # Create a vector containing the minimum and maximum values in data_matrix. 20 # Create a vector containing the minimum and maximum values in data_matrix.
30 min_max_vector = range(data_matrix); 21 min_max_vector = range(data_matrix);
31 # Create a color palette. 22 # Create a color palette.
85 } 76 }
86 state_color = get_state_color(data_matrix, histone_mark_color)[,2]; 77 state_color = get_state_color(data_matrix, histone_mark_color)[,2];
87 } 78 }
88 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); 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);
89 palette(defpalette); 80 palette(defpalette);
90 dev.off(); 81 if (!is.null(output_file_name)) {
82 dev.off();
83 }
84 return(state_color);
91 } 85 }
92 86
93 get_state_color <- function(data_matrix, histone_mark_color) { 87 get_state_color <- function(data_matrix, histone_mark_color) {
94 range_vector = apply(data_matrix, 1, range); 88 range_vector = apply(data_matrix, 1, range);
95 mm = NULL; 89 mm = NULL;
96 for(i in 1:dim(data_matrix)[1]) { 90 for(i in 1:dim(data_matrix)[1]) {
97 range_val1 = range_vector[1, i] + 1e-10 91 range_val1 = range_vector[1, i] + 1e-10;
98 range_val2 = range_vector[2, i] 92 range_val2 = range_vector[2, i];
99 mm = rbind(mm, (data_matrix[i,] - range_val1) / (range_val2 - range_val1)); 93 mm = rbind(mm, (data_matrix[i,] - range_val1) / (range_val2 - range_val1));
100 } 94 }
101 mm = mm^5; 95 mm = mm^5;
102 if(dim(mm)[2] > 1) { 96 if(dim(mm)[2] > 1) {
103 mm = mm / (apply(mm, 1, sum) + 1e-10); 97 mm = mm / (apply(mm, 1, sum) + 1e-10);
111 h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])})); 105 h = t(apply(state_color, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
112 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])}); 106 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
113 rt = cbind(rt, h); 107 rt = cbind(rt, h);
114 return(rt); 108 return(rt);
115 } 109 }
116
117 # Read the inputs.
118 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE);
119 for (i in 1:length(para_files)) {
120 para_file <- para_files[i];
121 para_file_base_name <- strsplit(para_file, split="/")[[1]][2]
122 output_file_base_name <- gsub(".para", "", para_file_base_name)
123 output_file_name <- paste(output_file_base_name, "state", i, "pdf", sep=".")
124 output_file_path <- paste(opt$output_dir, output_file_name, sep="/");
125 data_frame <- read.table(para_file, comment="!", header=T);
126 create_heatmap(data_frame, output_file_path);
127 }