annotate tools/mytools/intervalOverlap.py @ 3:33c067c3ae34

Deleted selected files
author xuebing
date Fri, 09 Mar 2012 19:47:53 -0500
parents 9071e359b9a3
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
1 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 find overlap and test signifiance
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 import os,sys
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7 def lineCount(filename):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 i = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 with open(filename) as f:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 for i, l in enumerate(f):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 return i + 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 def intersect(fileA,fileB,outfile,fraction,reciprocal):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 # return fileA intervals that overlap with interval in fileB
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 cmd = 'intersectBed -a '+fileA+' -b '+fileB + ' --wo -f '+fraction +' '+ reciprocal + '>'+outfile
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 #print cmd
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 os.system(cmd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 def parseIntersect(filename):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 # get number of overlapped A and B
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 nA = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 nB = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 return nA,nb
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 def shuffle(fileA,fileB,genomefile,fraction,reciprocal,N):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 # shuffle fileA N times, return the distribution of overlaps
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 nOverlap = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29 for i in range(N):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 # shuffle fileA using shuffleBed
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 #cmd = 'shuffleBed -i '+fileA+' -g '+genomefile +'>fileA.shuffled'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 # using random_interval.py
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33 cmd = 'python /Users/xuebing/galaxy-dist/tools/mytools/random_interval.py '+fileA+' fileA.shuffled across '+genomefile
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 os.system(cmd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 intersect('fileA.shuffled',fileB,'tmp',fraction,reciprocal)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 nOverlap.append(lineCount('tmp'))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 os.system('rm tmp')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 os.system('rm fileA.shuffled')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 return nOverlap
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 def main():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 fileA = sys.argv[1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 fileB = sys.argv[2]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 outfile = sys.argv[3]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 outplot = sys.argv[4]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 N = int(sys.argv[5]) # times to shuffle
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 genomefile = sys.argv[6]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 fraction = sys.argv[7]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 if len(sys.argv) == 9:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 reciprocal = sys.argv[8] # can only be '-r'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 reciprocal = ''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 print sys.argv
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 # intersect on real data
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 intersect(fileA,fileB,outfile,fraction,reciprocal)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 # number of overlaps
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 nOverlapReal = lineCount(outfile)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 print 'number of intervals in inputA that overlap with intervals in inputB:',nOverlapReal
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 # shuffle fileA to estimate background
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 nOverlapNull = shuffle(fileA,fileB,genomefile,fraction,reciprocal,N)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 # plot histogram
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 rscript = open('tmp.r','w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 rscript.write("x0 <- "+str(nOverlapReal)+"\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 rscript.write("x <- c("+','.join(map(str,nOverlapNull))+")\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 rscript.write("library(MASS)\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 rscript.write("\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 rscript.write("pv <- min((1+sum(x>=x0))/length(x),(1+sum(x<=x0))/length(x))\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 rscript.write("title <- paste('actual:chance = ',x0,':',round(mean(x)),' = ',format(x0/mean(x),digits=1,nsmall=2),', p-value < ',pv,sep='')\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 rscript.write("pdf('"+outplot+"')\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 rscript.write("h <- hist(x,breaks=50,xlab='number of overlaps',ylab='frequency',main=title)\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 rscript.write("plot(h$mids,h$counts,type='h',xlim=c(min(h$mids,x0),max(x0,h$mids)),ylim=c(0,max(h$counts)),xlab='number of overlaps',ylab='frequency',main=title)\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 rscript.write("points(x0,0,col='red')\n")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 rscript.write("dev.off()")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 rscript.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 os.system("R --vanilla < tmp.r")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 os.system('rm tmp.r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 main()