|
105
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
|
4
|
|
|
5 option_list <- list(
|
|
121
|
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")
|
|
105
|
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
|
|
133
|
14 create_heatmap<-function(data_frame, output_file_name) {
|
|
|
15 # Plot a heatmap for a .para / .state combination
|
|
|
16 # based on the received data_frame which was created
|
|
|
17 # by reading the .para file.
|
|
|
18 num_columns = dim(data_frame)[2];
|
|
|
19 num_rows = dim(data_frame)[1];
|
|
|
20 p = (sqrt(9 + 8 * (num_columns - 1)) - 3) / 2;
|
|
|
21 data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
|
|
|
22 column_names(data_matrix) = column_names(data_frame)[1+1:p];
|
|
|
23 marks = column_names(data_matrix);
|
|
|
24 row_names(data_matrix) = paste(1:num_rows, " (", round(data_frame[,1]/sum(data_frame[,1])*10000)/100, "%)", sep="");
|
|
|
25
|
|
|
26 # Open the output PDF file.
|
|
|
27 pdf(file=output_file_name);
|
|
|
28 # Set graphical parameters.
|
|
|
29 par(mar=c(6, 1, 1, 6));
|
|
|
30 # Create a vector containing the minimum and maximum values in data_matrix.
|
|
|
31 min_max_vector = range(data_matrix);
|
|
|
32 # Define colors for the palette.
|
|
|
33 colors = 0:100 / 100 * (min_max_vector[2] - min_max_vector[1]) + min_max_vector[1];
|
|
|
34 # Create the color palette.
|
|
|
35 my_palette = colorRampPalette(c("white", "dark blue"))(n=100);
|
|
|
36 defpalette = palette(my_palette);
|
|
|
37 # Plot the heatmap for the current .para / .state combination.
|
|
|
38 plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
|
|
|
39 axis(1, at=1:p-0.5, labels=column_names(data_matrix), las=2);
|
|
|
40 axis(4, at=1:num_rows-0.5, labels=row_names(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),
|
|
|
42 col=round((t(data_matrix)-min_max_vector[1])/(min_max_vector[2]-min_max_vector[1])*100));
|
|
|
43 markcolor = t(col2rgb(terrain.colors(ceiling(p))[1:p]));
|
|
|
44
|
|
|
45 for(i in 1:length(marks)) {
|
|
|
46 if(regexpr("h3k4me3", tolower(marks[i])) > 0) {
|
|
|
47 markcolor[i,] = c(255, 0, 0);
|
|
|
48 }
|
|
|
49 if(regexpr("h3k4me2", tolower(marks[i])) > 0) {
|
|
|
50 markcolor[i,] = c(250, 100, 0);
|
|
|
51 }
|
|
|
52 if(regexpr("h3k4me1", tolower(marks[i])) > 0) {
|
|
|
53 markcolor[i,] = c(250, 250, 0);
|
|
|
54 }
|
|
|
55 if(regexpr("h3k36me3", tolower(marks[i]))>0) {
|
|
|
56 markcolor[i,] = c(0, 150, 0);
|
|
|
57 }
|
|
|
58 if(regexpr("h2a", tolower(marks[i])) > 0) {
|
|
|
59 markcolor[i,] = c(0, 150, 150);
|
|
|
60 }
|
|
|
61 if(regexpr("dnase", tolower(marks[i])) > 0) {
|
|
|
62 markcolor[i,] = c(0, 200, 200);
|
|
|
63 }
|
|
|
64 if(regexpr("h3k9ac", tolower(marks[i])) > 0) {
|
|
|
65 markcolor[i,] = c(250, 0, 200);
|
|
|
66 }
|
|
|
67 if(regexpr("h3k9me3", tolower(marks[i])) > 0) {
|
|
|
68 markcolor[i,] = c(100, 100, 100);
|
|
|
69 }
|
|
|
70 if(regexpr("h3k27ac", tolower(marks[i])) > 0) {
|
|
|
71 markcolor[i,] = c(250, 150, 0);
|
|
|
72 }
|
|
|
73 if(regexpr("h3k27me3", tolower(marks[i])) > 0) {
|
|
|
74 markcolor[i,] = c(0, 0, 200);
|
|
|
75 }
|
|
|
76 if(regexpr("h3k79me2", tolower(marks[i])) > 0) {
|
|
|
77 markcolor[i,] = c(200, 0, 200);
|
|
|
78 }
|
|
|
79 if(regexpr("h4k20me1", tolower(marks[i])) > 0) {
|
|
|
80 markcolor[i,] = c(50, 200, 50);
|
|
|
81 }
|
|
|
82 if(regexpr("ctcf", tolower(marks[i])) > 0) {
|
|
|
83 markcolor[i,] = c(200, 0, 250);
|
|
|
84 }
|
|
|
85 state_color = get_state_color(data_matrix, markcolor)[,2];
|
|
118
|
86 }
|
|
133
|
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);
|
|
|
88 palette(defpalette);
|
|
|
89 dev.off();
|
|
|
90 }
|
|
118
|
91
|
|
133
|
92 get_state_color<-function(statemean, markcolor) {
|
|
|
93 rg=apply(statemean, 1, range);
|
|
|
94 mm=NULL;
|
|
118
|
95 for(i in 1:dim(statemean)[1]) {
|
|
133
|
96 mm = rbind(mm, (statemean[i,]-rg[1,i]+1e-10)/(rg[2,i]-rg[1,i]+1e-10));
|
|
118
|
97 }
|
|
133
|
98 mm = mm^5;
|
|
118
|
99 if(dim(mm)[2] > 1) {
|
|
|
100 mm = mm / (apply(mm, 1, sum) + 1e-10);
|
|
|
101 }
|
|
|
102 mycol = mm%*%markcolor;
|
|
|
103 s = apply(statemean, 1, max);
|
|
133
|
104 s = (s-min(s)) / (max(s)-min(s) + 1e-10);
|
|
|
105 mycol = round(255-(255-mycol)*s/0.5);
|
|
|
106 mycol[mycol<0] = 0;
|
|
|
107 rt = paste(mycol[,1], mycol[,2], mycol[,3], sep=",");
|
|
118
|
108 h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
|
|
|
109 h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
|
|
133
|
110 rt = cbind(rt, h);
|
|
118
|
111 return(rt);
|
|
|
112 }
|
|
|
113
|
|
121
|
114 # Read the inputs.
|
|
122
|
115 para_files <- list.files(path=opt$input_dir, pattern="\\.para$", full.names=TRUE);
|
|
|
116 for (i in 1:length(para_files)) {
|
|
121
|
117 para_file <- para_files[i];
|
|
122
|
118 para_file_base_name <- strsplit(para_file, split="/")[[1]][2]
|
|
|
119 output_file_name <- gsub(".para", ".pdf", para_file_base_name)
|
|
133
|
120 output_file_path <- paste(opt$output_dir, output_file_name, sep="/");
|
|
121
|
121 data_frame <- read.table(para_file, comment="!", header=T);
|
|
122
|
122 create_heatmap(data_frame, output_file_path);
|
|
121
|
123 }
|
