annotate mytools/alignvis.py @ 0:39217fa39ff2

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