|
105
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3 suppressPackageStartupMessages(library("optparse"))
|
|
|
4
|
|
|
5 option_list <- list(
|
|
|
6 make_option(c("-i", "--input"), action="store", dest="input", help="IDEAS para file"),
|
|
|
7 make_option(c("-o", "--output"), action="store", dest="output", help="Output PDF file")
|
|
|
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_matrix, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) {
|
|
|
15 k = dim(data_matrix)[2];
|
|
|
16 l = dim(data_matrix)[1];
|
|
|
17 p = (sqrt(9 + 8 * (k - 1)) - 3) / 2;
|
|
|
18 m = as.matrix(data_matrix[,1+1:p]/data_matrix[,1]);
|
|
114
|
19 colnames(m) = colnames(data_matrix)[1+1:p];
|
|
105
|
20 marks = colnames(m);
|
|
|
21 rownames(m) = paste(1:l-1," (", round(data_matrix[,1] / sum(data_matrix[,1]) * 10000) / 100, "%)", sep="");
|
|
|
22 pdf(file=opt$output);
|
|
|
23 par(mar=c(6,1,1,6));
|
|
|
24 rg = range(m);
|
|
|
25 colors = 0:100/100*(rg[2]-rg[1])+rg[1];
|
|
|
26 my_palette = colorRampPalette(cols)(n=100);
|
|
|
27 defpalette = palette(my_palette);
|
|
|
28
|
|
|
29 plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
|
|
|
30 axis(1, at=1:p-0.5, labels=colnames(m), las=2);
|
|
|
31 axis(4, at=1:l-0.5, labels=rownames(m), las=2);
|
|
|
32 rect(rep(1:p-1,l), rep(1:l-1,each=p), rep(1:p,l), rep(1:l,each=p), col=round((t(m)-rg[1])/(rg[2]-rg[1])*100));
|
|
|
33
|
|
|
34 if(length(statecolor)==0) {
|
|
|
35 if(length(markcolor)==0) {
|
|
|
36 markcolor=t(col2rgb(terrain.colors(ceiling(p))[1:p]));
|
|
|
37 for(i in 1:length(marks)) {
|
|
|
38 if(regexpr("h3k4me3",tolower(marks[i]))>0) {
|
|
|
39 markcolor[i,]=c(255,0,0);
|
|
|
40 }
|
|
|
41 if(regexpr("h3k4me2",tolower(marks[i]))>0) {
|
|
|
42 markcolor[i,]=c(250,100,0);
|
|
|
43 }
|
|
|
44 if(regexpr("h3k4me1",tolower(marks[i]))>0) {
|
|
|
45 markcolor[i,]=c(250,250,0);
|
|
|
46 }
|
|
|
47 if(regexpr("h3k36me3",tolower(marks[i]))>0) {
|
|
|
48 markcolor[i,]=c(0,150,0);
|
|
|
49 }
|
|
|
50 if(regexpr("h2a",tolower(marks[i]))>0) {
|
|
|
51 markcolor[i,]=c(0,150,150);
|
|
|
52 }
|
|
|
53 if(regexpr("dnase",tolower(marks[i]))>0) {
|
|
|
54 markcolor[i,]=c(0,200,200);
|
|
|
55 }
|
|
|
56 if(regexpr("h3k9ac",tolower(marks[i]))>0) {
|
|
|
57 markcolor[i,]=c(250,0,200);
|
|
|
58 }
|
|
|
59 if(regexpr("h3k9me3",tolower(marks[i]))>0) {
|
|
|
60 markcolor[i,]=c(100,100,100);
|
|
|
61 }
|
|
|
62 if(regexpr("h3k27ac",tolower(marks[i]))>0) {
|
|
|
63 markcolor[i,]=c(250,150,0);
|
|
|
64 }
|
|
|
65 if(regexpr("h3k27me3",tolower(marks[i]))>0) {
|
|
|
66 markcolor[i,]=c(0,0,200);
|
|
|
67 }
|
|
|
68 if(regexpr("h3k79me2",tolower(marks[i]))>0) {
|
|
|
69 markcolor[i,]=c(200,0,200);
|
|
|
70 }
|
|
|
71 if(regexpr("h4k20me1",tolower(marks[i]))>0) {
|
|
|
72 markcolor[i,]=c(50,200,50);
|
|
|
73 }
|
|
|
74 if(regexpr("ctcf",tolower(marks[i]))>0) {
|
|
|
75 markcolor[i,]=c(200,0,250);
|
|
|
76 }
|
|
|
77 }
|
|
|
78 }
|
|
|
79 statecolor = stateColor(m, markcolor)[,2];
|
|
|
80 }
|
|
|
81 rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=statecolor);
|
|
|
82
|
|
|
83 palette(defpalette);
|
|
|
84 dev.off();
|
|
|
85 #return(statecolor);
|
|
|
86 }
|
|
|
87
|
|
|
88 # Read the input.
|
|
|
89 data_matrix <- read.table(opt$input, header=T);
|
|
|
90 create_heatmap(data_matrix);
|
