[3] | 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() |
---|