annotate getGenomicScore.py @ 12:2f4ea569f048

Uploaded
author xuebing
date Sat, 10 Mar 2012 08:10:44 -0500
parents b7f1d9f8f3bc
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
2 import random,string,os,sys
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
3
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
4
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
5 def getScore(intvfile,outfile,summary_type,bwfilepath,nbin,strand,outplot,span):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
6 f = open(intvfile)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
7 tmpsh = "".join(random.sample(string.letters+string.digits, 8))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
8 tmpout = "".join(random.sample(string.letters+string.digits, 8))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
9 tmp = open(tmpsh,'w')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
10 if os.path.isdir(bwfilepath):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
11 for line in f:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
12 flds = line.strip().split('\t')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
13 cmd = 'bigWigSummary -type='+summary_type+' '+bwfilepath+'/'+flds[0]+'.bw '+flds[0]+' '+flds[1]+' '+flds[2]+' '+nbin+' >> '+tmpout+' 2>>'+tmpout+'\n'
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
14 tmp.write(cmd)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
15 else:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
16 for line in f:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
17 flds = line.strip().split('\t')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
18 cmd = 'bigWigSummary -type='+summary_type+' '+bwfilepath+' '+flds[0]+' '+flds[1]+' '+flds[2]+' '+nbin+' >> '+tmpout+' 2>>'+tmpout+'\n'
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
19 tmp.write(cmd)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
20 f.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
21 # remove blank lines
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
22 tmp.write("sed '/^$/d' "+tmpout+'>'+tmpout+".1\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
23 tmp.write("sed '/^Can/d' "+tmpout+".1 >"+tmpout+".2\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
24 # set n/a to 0
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
25 tmp.write("sed 's/n\/a/0/g' "+tmpout+".2 >"+tmpout+".3\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
26 # replace text with 0
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
27 zeros = ''.join(['0\t']*int(nbin))
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
28 tmp.write("sed 's/^[a-zA-Z]/"+zeros+"/' "+tmpout+".3 >"+tmpout+".4\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
29 # cut the first nbin columns
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
30 tmp.write("cut -f 1-"+nbin+" "+tmpout+".4 > "+tmpout+".5\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
31 tmp.write("paste "+intvfile+" "+tmpout+".5 >"+outfile+"\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
32 tmp.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
33 os.system('chmod +x '+tmpsh)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
34 os.system('./'+tmpsh)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
35 #os.system('rm '+tmpout+'*')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
36 #os.system('rm '+tmpsh)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
37
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
38 # strandness: need to reverse bins for - strand
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
39 if nbin > 1 and strand > 0:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
40 strand = strand - 1 # python is 0 based
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
41 os.system('mv '+outfile+' '+tmpout)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
42 f = open(tmpout)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
43 out = open(outfile,'w')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
44 for line in f:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
45 flds=line.strip().split('\t')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
46 if flds[strand] == '+':
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
47 out.write(line)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
48 else:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
49 scores = flds[-int(nbin):]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
50 scores.reverse()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
51 flds = flds[:-int(nbin)]+scores
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
52 out.write('\t'.join(flds)+'\n')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
53 os.system('rm '+tmpout)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
54 f.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
55 out.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
56 # plot
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
57 if int(nbin) > 1:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
58 rscript = open(tmpsh,"w")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
59 rscript.write("options(warn=-1)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
60 rscript.write("x <- read.table('"+outfile+"',sep='\t')\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
61 rscript.write("x <- x[,(ncol(x)+1-"+nbin+"):ncol(x)]\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
62 rscript.write("pdf('"+outplot+"')\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
63 rscript.write("avg <- apply(x,2,mean)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
64 rscript.write("err <- apply(x,2,sd)/sqrt(nrow(x))\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
65 rscript.write("ylim=c(min(avg-err),max(avg+err))\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
66 rscript.write("xticks <- seq(ncol(x))-(1+ncol(x))/2\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
67 if span >= 0.1:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
68 rscript.write("avg = loess(avg~xticks,span="+str(span)+")$fitted\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
69 rscript.write("err = loess(err~xticks,span="+str(span)+")$fitted\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
70 rscript.write("par(cex=1.5)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
71 rscript.write("plot(xticks,avg,ylab='average conservation score',xlab='relative position (bin)',type='l',lwd=0,ylim=ylim)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
72 rscript.write("polygon(c(xticks,rev(xticks)),c(avg+err,rev(avg-err)),col='slateblue1',border=NA)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
73 rscript.write("lines(xticks,avg,type='l',lwd=1)\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
74 rscript.write("dev.off()\n")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
75 rscript.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
76 os.system("R --vanilla < "+tmpsh)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
77 os.system("rm "+tmpsh)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
78
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
79 getScore(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],int(sys.argv[6]),sys.argv[7],float(sys.argv[8]))