Mercurial > repos > xuebing > sharplab_seq_motif
view mytools/alignvis.py @ 0:39217fa39ff2
Uploaded
author | xuebing |
---|---|
date | Tue, 13 Mar 2012 23:34:52 -0400 |
parents | |
children |
line wrap: on
line source
import sys,os infile = sys.argv[1] outfile = sys.argv[2] uselog = sys.argv[3] subset = sys.argv[4] reorder = sys.argv[5] color = sys.argv[6] scale = sys.argv[7] # rescale each row rscript = open('tmp.r','w') rscript.write("x <- read.table('"+infile+"')\n") rscript.write("nfeat <- nrow(x) \n") rscript.write("nbin <- ncol(x) - 2\n") rscript.write("totalcov <- x[,2]\n") rscript.write("x <- x[,3:ncol(x)]\n") if subset =='subset': rscript.write("if (nfeat*nbin > 100000) {\n") rscript.write(" nfeat2 <- as.integer(100000/nbin)\n") rscript.write(" subind <- sample(seq(nfeat),nfeat2)\n") rscript.write(" x <- x[subind,]\n") rscript.write(" totalcov <- totalcov[subind]\n") rscript.write("}\n") rscript.write("pdf('"+outfile+"')\n") if uselog == 'uselog': rscript.write("x <- -(log(1+as.matrix(x,nc=ncol(x)-2)))\n") else: rscript.write("x <- -as.matrix(x,nc=ncol(x)-2)\n") if scale == 'scale': rscript.write("x <- scale(x)\n") if reorder == 'average': rscript.write("hc <- hclust(dist(x),method= 'average')\n") rscript.write("x <- x[hc$order,]\n") elif reorder == 'centroid': rscript.write("hc <- hclust(dist(x),method= 'centroid')\n") rscript.write("x <- x[hc$order,]\n") elif reorder == 'complete': rscript.write("hc <- hclust(dist(x),method= 'complete')\n") rscript.write("x <- x[hc$order,]\n") elif reorder == 'single': rscript.write("hc <- hclust(dist(x),method= 'single')\n") rscript.write("x <- x[hc$order,]\n") elif reorder == 'median': rscript.write("hc <- hclust(dist(x),method= 'median')\n") rscript.write("x <- x[hc$order,]\n") elif reorder == 'sort_by_total': rscript.write("srt <- sort(totalcov,index.return=T)\n") rscript.write("x <- x[srt$ix,]\n") elif reorder == 'sort_by_center': rscript.write("srt <- sort(x[,as.integer(nbin/2)],index.return=T)\n") rscript.write("x <- x[srt$ix,]\n") if color == 'heat': rscript.write("colormap = heat.colors(1000)\n") elif color == 'topo': rscript.write("colormap = topo.colors(1000)\n") elif color == 'rainbow': rscript.write("colormap = rainbow(1000)\n") elif color == 'terrain': rscript.write("colormap = terrain.colors(1000)\n") else: rscript.write("colormap = gray.colors(1000)\n") #rscript.write("qt <- quantile(as.vector(x),probs=c(0.1,0.9))\n") #rscript.write("breaks <- c(min(as.vector(x)),seq(qt[1],qt[2],length.out=99),max(as.vector(x)))\n") #rscript.write("image(t(x),col=colormap,breaks=breaks,axes=F)\n") rscript.write("image(t(x),col=colormap,axes=F)\n") rscript.write("dev.off()\n") rscript.close() os.system("R --slave < tmp.r") os.system("rm tmp.r")