root/galaxy-central/eggs/bx_python-0.5.0_dev_f74aec067563-py2.6-macosx-10.6-universal-ucs2.egg/bx/pwm/maf_select_motifs.py

#!/usr/bin/env python2.4
"""
Returns all positions of a maf with any pwm score > threshold
The positions are projected onto human coordinates
"""

import psyco_full
from bx.align import maf as align_maf
import bx.pwm.position_weight_matrix as pwmx
from bx.pwm.pwm_score_maf import MafMotifSelect
import sys
from bx import intervals

def isnan(x):
    return not x==x

def main():

    if len(sys.argv) < 5:
        print >>sys.stderr, "%s transfac|basic pwmfile inmaf threshold [motif]" % sys.argv[0]
        sys.exit(2)

    r = pwmx.Reader(open(sys.argv[2]),format=sys.argv[1])
    pwm = iter(r).next()
    inmaf = open(sys.argv[3])
    threshold = float(sys.argv[4])
    if len(sys.argv) > 5: motif = sys.argv[5]
    else: motif = None

    for maf in align_maf.Reader(inmaf):
        for mafmotif,pwm_score,motif_score in MafMotifSelect(maf, pwm, motif, threshold):
            #mafwrite( mafmotif,pwm_score,motif_score)
            print mafmotif, pwm_score, motif_score
            print 'zzzzzzzzzzzzzzzzzzzzzzzzzzzzz'

def mafwrite(alignment,kvec=None,jvec=None,file=sys.stdout):
    file.write( "a score=" + str( alignment.score ) )
    for key in alignment.attributes:
        file.write( " %s=%s" % ( key, alignment.attributes[key] ) )
    file.write( "\n" )
    rows = []
    if not kvec: kvec = [0 for c in alignment.components]
    if not jvec: jvec = [0 for c in alignment.components]
    for c,x,y in zip(alignment.components,kvec,jvec):
    #for c in alignment.components:
        rows.append( ( "s", c.src, str( c.start ), str( c.size ), c.strand, str( c.src_size ), c.text, "%.2f" % x, str(y) ) )
        #rows.append( ( "s", c.src, str( c.start ), str( c.size ), c.strand, str( c.src_size ), c.text ) )
        file.write( format_tabular( rows, "llrrrrrrr" ) )

def format_tabular( rows, align=None ):
    if len( rows ) == 0: return ""
    lengths = [ len( col ) for col in rows[ 0 ] ]
    for row in rows[1:]:
        for i in range( 0, len( row ) ):
            lengths[ i ] = max( lengths[ i ], len( row[ i ] ) )
    rval = ""
    for row in rows:
        for i in range( 0, len( row ) ):
            if align and align[ i ] == "l":
                rval += row[ i ].ljust( lengths[ i ] )
            else:
                rval += row[ i ].rjust( lengths[ i ] )
            rval += " "
        rval += "\n"
    return rval

    

if __name__ == '__main__': main()