|
11
|
1 import sys,os
|
|
|
2
|
|
|
3 infile = sys.argv[1]
|
|
|
4 outfile = sys.argv[2]
|
|
|
5 uselog = sys.argv[3]
|
|
|
6 subset = sys.argv[4]
|
|
|
7 reorder = sys.argv[5]
|
|
|
8 color = sys.argv[6]
|
|
|
9 scale = sys.argv[7] # rescale each row
|
|
|
10 cols = sys.argv[8]
|
|
|
11 rscript = open('tmp.r','w')
|
|
|
12
|
|
|
13 rscript.write("x <- read.table('"+infile+"')\n")
|
|
|
14 rscript.write("x <- x[,c("+cols+")]\n")
|
|
|
15 rscript.write("nr <- nrow(x) \n")
|
|
|
16 rscript.write("nc <- ncol(x)\n")
|
|
|
17 rscript.write("rowsum <- apply(x,1,sum)\n")
|
|
|
18
|
|
|
19 if subset =='subset':
|
|
|
20 rscript.write("if (nr*nc > 100000) {\n")
|
|
|
21 rscript.write(" nr2 <- as.integer(100000/nc)\n")
|
|
|
22 rscript.write(" subind <- sample(seq(nr),nr2)\n")
|
|
|
23 rscript.write(" x <- x[subind,]\n")
|
|
|
24 rscript.write(" rowsum <- rowsum[subind]\n")
|
|
|
25 rscript.write("}\n")
|
|
|
26
|
|
|
27 rscript.write("pdf('"+outfile+"')\n")
|
|
|
28
|
|
|
29 if uselog == 'uselog':
|
|
|
30 rscript.write("x <- -(log(as.matrix(x,nc=nc)))\n")
|
|
|
31 else:
|
|
|
32 rscript.write("x <- -as.matrix(x,nc=nc)\n")
|
|
|
33 if scale == 'scale':
|
|
|
34 rscript.write("x <- scale(x)\n")
|
|
|
35 if reorder == 'average':
|
|
|
36 rscript.write("hc <- hclust(dist(x),method= 'average')\n")
|
|
|
37 rscript.write("x <- x[hc$order,]\n")
|
|
|
38 elif reorder == 'centroid':
|
|
|
39 rscript.write("hc <- hclust(dist(x),method= 'centroid')\n")
|
|
|
40 rscript.write("x <- x[hc$order,]\n")
|
|
|
41 elif reorder == 'complete':
|
|
|
42 rscript.write("hc <- hclust(dist(x),method= 'complete')\n")
|
|
|
43 rscript.write("x <- x[hc$order,]\n")
|
|
|
44 elif reorder == 'single':
|
|
|
45 rscript.write("hc <- hclust(dist(x),method= 'single')\n")
|
|
|
46 rscript.write("x <- x[hc$order,]\n")
|
|
|
47 elif reorder == 'median':
|
|
|
48 rscript.write("hc <- hclust(dist(x),method= 'median')\n")
|
|
|
49 rscript.write("x <- x[hc$order,]\n")
|
|
|
50 elif reorder == 'sort_by_total':
|
|
|
51 rscript.write("srt <- sort(rowsum,index.return=T)\n")
|
|
|
52 rscript.write("x <- x[srt$ix,]\n")
|
|
|
53 elif reorder == 'sort_by_center':
|
|
|
54 rscript.write("srt <- sort(x[,as.integer(nc/2)],index.return=T)\n")
|
|
|
55 rscript.write("x <- x[srt$ix,]\n")
|
|
|
56 if color == 'heat':
|
|
|
57 rscript.write("colormap = heat.colors(1000)\n")
|
|
|
58 elif color == 'topo':
|
|
|
59 rscript.write("colormap = topo.colors(1000)\n")
|
|
|
60 elif color == 'rainbow':
|
|
|
61 rscript.write("colormap = rainbow(1000)\n")
|
|
|
62 elif color == 'terrain':
|
|
|
63 rscript.write("colormap = terrain.colors(1000)\n")
|
|
|
64 else:
|
|
|
65 rscript.write("colormap = gray.colors(1000)\n")
|
|
|
66
|
|
|
67 #rscript.write("qt <- quantile(as.vector(x),probs=c(0.1,0.9))\n")
|
|
|
68 #rscript.write("breaks <- c(min(as.vector(x)),seq(qt[1],qt[2],length.out=99),max(as.vector(x)))\n")
|
|
|
69 #rscript.write("image(t(x),col=colormap,breaks=breaks,axes=F)\n")
|
|
|
70 rscript.write("image(t(x),col=colormap,axes=F)\n")
|
|
|
71 rscript.write("dev.off()\n")
|
|
|
72
|
|
|
73 rscript.close()
|
|
|
74
|
|
|
75 os.system("R --slave < tmp.r")
|
|
|
76 os.system("rm tmp.r")
|
|
|
77
|
