Mercurial > repos > xuebing > sharplabtool
diff align2database.py @ 14:76e1b1b21cce default tip
Deleted selected files
author | xuebing |
---|---|
date | Tue, 13 Mar 2012 19:05:10 -0400 |
parents | 292186c14b08 |
children |
line wrap: on
line diff
--- a/align2database.py Sat Mar 10 08:17:36 2012 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,104 +0,0 @@ -''' -align mulitple bed to one bed -python align_multiple.py hmChIP_mm9_peak_bed/mm9-GSE19999_PolII_P25_all.cod.bed hmChIP_mm9_peak_bed/ test.txt test.pdf 100 5000 -''' - -import os,sys,random -def main(): - queryfile = sys.argv[1] - inpath = sys.argv[2] - outputdata = sys.argv[3] - outputerr = sys.argv[4] - barplotpdf = sys.argv[5] - min_feat = sys.argv[6] # min features overlap - windowsize = sys.argv[7] - anchor = sys.argv[8] - span = sys.argv[9] # loess smooth parameter - - inpath = inpath.rstrip('/') - #os.system('rm '+inpath+'/*tmp*') - - infiles = os.listdir(inpath) - - #print len(infiles),' files\n' - i = 0 - for infile in infiles: - if 'tmp' in infile: - #os.system('rm '+inpath+'/'+infile) - continue - i = i +1 - print i,infile - output = infile.split('/')[-1]+'-to-'+queryfile.split('/')[-1]#'.summary' - if anchor == 'database': - command = 'python /Users/xuebing/galaxy-dist/tools/mytools/alignr.py -b '+inpath+'/'+infile+' -a '+queryfile+' -o '+output+' --summary-only -q -w '+windowsize - else: - command = 'python /Users/xuebing/galaxy-dist/tools/mytools/alignr.py -a '+inpath+'/'+infile+' -b '+queryfile+' -o '+output+' --summary-only -q -w '+windowsize - #print command+'\n' - os.system(command) - print 'start visualization...' - # visualize - rscriptfile = 'f-'+str(random.random())+'.r' - r = open(rscriptfile,'w') - r.write("files <- dir('.','summary',full.name=T)\n") - #r.write("print(files)\n") - r.write("x <- read.table(files[1])\n") - r.write("err <- read.table(gsub('summary','standarderror',files[1]))\n") - r.write("for (filename in files[2:length(files)]){\n") - r.write(" x <- rbind(x,read.table(filename))\n") - r.write(" err <- rbind(err,read.table(gsub('summary','standarderror',filename)))\n") - r.write("}\n") - r.write("x <- x[x[,2] > "+min_feat+",]\n") - r.write("err <- err[err[,2] > "+min_feat+",]\n") - r.write("name <- as.character(x[,1])\n") - r.write("nfeat <- x[,2]\n") - r.write("x <- x[,3:ncol(x)]\n") - r.write("err <- err[,3:ncol(err)]\n") - r.write("for (i in 1:nrow(x)) {\n") - r.write(" name[i] <- strsplit(name[i],'-to-')[[1]][1]\n") - r.write("}\n") - # remove rows that have no variation, which cause problem in heatmap. This is the case when align to itself - r.write("toremove <- seq(nrow(x))\n") - r.write("for (i in 1:nrow(x)){\n") - r.write(" toremove[i] <- all(x[i,] == x[i,1])\n") - r.write("}\n") - r.write("x <- x[!toremove,]\n") - r.write("err <- err[!toremove,]\n") - r.write("name <- name[!toremove]\n") - r.write("nfeat <- nfeat[!toremove]\n") - r.write("write.table(cbind(name,nfeat,x),file='"+outputdata+"',sep ='\\t',quote=F,row.names=F,col.names=F)\n") - r.write("write.table(cbind(name,nfeat,err),file='"+outputerr+"',sep ='\\t',quote=F,row.names=F,col.names=F)\n") - - r.write("pdf('"+barplotpdf+"')\n") - r.write("heatmap(t(scale(t(as.matrix(x,nc=ncol(x))))),Colv=NA,labRow=name,labCol=NA,margins=c(1,8),cexRow=0.02+1/log(nrow(x)),col=heat.colors(1000))\n") - - if windowsize != '0' : - r.write("xticks <- seq(-"+windowsize+","+windowsize+",length.out=100)\n") - r.write("xlab <- 'relative position (bp)'\n") - else: - r.write("xticks <- seq(100)\n") - r.write("xlab <- 'relative position (bin)'\n") - - #r.write("plot(xticks,colSums(t(scale(t(as.matrix(x,nc=ncol(x)))))),xlab='relative position (bp)',type='l',lwd=2,main='total signal')\n") - r.write("for (i in 1:nrow(x)) {\n") - r.write(" avg <- x[i,]/nfeat[i]\n") - #r.write(" maxv <- max(avg)\n") - #r.write(" minv <- min(avg)\n") - #r.write(" medv <- median(avg)\n") - #r.write(" if (maxv < "+fold+"*medv | minv*"+fold+">medv){next}\n") - #smooth - if float(span) >= 0.1: - r.write(" avg = loess(as.numeric(avg)~xticks,span="+span+")$fitted\n") - r.write(" err[i,] = loess(as.numeric(err[i,])~xticks,span="+span+")$fitted\n") - r.write(" par(cex=1.5)\n") - r.write(" plot(xticks,avg,ylab='average coverage',main=paste(name[i],'\n n=',nfeat[i],sep=''),xlab=xlab,type='l',lwd=1,ylim=c(min(avg-err[i,]),max(avg+err[i,])))\n") - r.write(" polygon(c(xticks,rev(xticks)),c(avg+err[i,],rev(avg-err[i,])),col='slateblue1',border=NA)\n") - r.write(" lines(xticks,avg,type='l',lwd=1)\n") - r.write("}\n") - r.write("dev.off()\n") - r.close() - os.system("R --vanilla < "+rscriptfile) - os.system('rm '+rscriptfile) - os.system('rm *.summary') - os.system('rm *.standarderror') - -main()