comparison mytools/alignvis.py @ 0:39217fa39ff2

Uploaded

0:39217fa39ff2

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")