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