1 | #!/usr/bin/env python2.4 |
---|
2 | """ |
---|
3 | Returns all positions of a maf with any pwm score > threshold |
---|
4 | The positions are projected onto human coordinates |
---|
5 | """ |
---|
6 | |
---|
7 | import psyco_full |
---|
8 | from bx.align import maf as align_maf |
---|
9 | import position_weight_matrix as pwmx |
---|
10 | from bx.pwm.pwm_score_maf import MafMotifScorer |
---|
11 | import sys |
---|
12 | from bx import intervals |
---|
13 | |
---|
14 | def isnan(x): |
---|
15 | return not x==x |
---|
16 | |
---|
17 | def main(): |
---|
18 | |
---|
19 | if len(sys.argv) < 4: |
---|
20 | print >>sys.stderr, "%s motif inmaf spec1,spec2,... " % sys.argv[0] |
---|
21 | sys.exit(0) |
---|
22 | |
---|
23 | targmotif = sys.argv[1] |
---|
24 | inmaf = open(sys.argv[2]) |
---|
25 | threshold = 0 |
---|
26 | |
---|
27 | species = [] |
---|
28 | |
---|
29 | for sp in sys.argv[3].split(','): |
---|
30 | species.append( sp ) |
---|
31 | |
---|
32 | for maf in align_maf.Reader(inmaf): |
---|
33 | mafchrom = maf.components[0].src.split('.')[1] |
---|
34 | mafstart = maf.components[0].start |
---|
35 | mafend = maf.components[0].end |
---|
36 | reftext = maf.components[0].text |
---|
37 | |
---|
38 | # maf block scores for each matrix |
---|
39 | for scoremax,width,headers in MafMotifScorer(species, maf,targmotif): |
---|
40 | #print >>sys.stderr,headers |
---|
41 | blocklength = width |
---|
42 | mafsrc,mafstart,mafend = headers[0] |
---|
43 | mafchrom = mafsrc.split('.')[1] |
---|
44 | |
---|
45 | # lists of scores for each position in scoremax |
---|
46 | mx = scoremax |
---|
47 | for offset in range(blocklength): |
---|
48 | |
---|
49 | # scan all species with threshold |
---|
50 | for i in range(len(species)): |
---|
51 | if mx[i][offset] > threshold: |
---|
52 | refstart = mafstart + offset - reftext.count('-',0,offset) |
---|
53 | refend = refstart + len( targmotif ) |
---|
54 | data = " ".join([ "%.2f" % mx[x][offset] for x in range(len(species))]) |
---|
55 | # quote the motif |
---|
56 | print mafchrom,refstart,refend,"'"+targmotif+"'",data |
---|
57 | break |
---|
58 | |
---|
59 | if __name__ == '__main__': main() |
---|