/tools/meme/fimo_wrapper.py
https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 73 lines · 54 code · 9 blank · 10 comment · 13 complexity · 53354ffd0f8dcd474aa4455b7fefd8e2 MD5 · raw file
- #!/usr/bin/env python
- #Dan Blankenberg
- """
- Read text output from FIMO and create an interval file.
- """
- import sys, tempfile, subprocess, shutil, os
- from galaxy_utils.sequence.transform import DNA_reverse_complement
- buffsize = 1048576
- def stop_err( msg ):
- sys.stderr.write( msg )
- sys.exit()
- def main():
- assert len( sys.argv ) == 8, "Wrong number of arguments"
- sys.argv.pop(0)
- fimo_cmd = sys.argv.pop(0)
- html_path = sys.argv.pop(0)
- html_out = sys.argv.pop(0)
- interval_out = sys.argv.pop(0)
- txt_out = sys.argv.pop(0)
- xml_out = sys.argv.pop(0)
- gff_out = sys.argv.pop(0)
-
- #run fimo
- try:
- tmp_stderr = tempfile.NamedTemporaryFile()
- #tmp_stderr = open( tmp_filename, 'wb' )
- proc = subprocess.Popen( args=fimo_cmd, shell=True, stderr=tmp_stderr )
- returncode = proc.wait()
- #tmp_stderr.close()
- # get stderr, allowing for case where it's very large
- #tmp_stderr = open( tmp, 'rb' )
- tmp_stderr.seek(0)
- stderr = ''
- try:
- while True:
- stderr += tmp_stderr.read( buffsize )
- if not stderr or len( stderr ) % buffsize != 0:
- break
- except OverflowError:
- pass
-
- if returncode != 0:
- raise Exception, stderr
- except Exception, e:
- raise Exception, 'Error running FIMO:\n' + str( e )
- shutil.move( os.path.join( html_path, 'fimo.txt' ), txt_out )
- shutil.move( os.path.join( html_path, 'fimo.gff' ), gff_out )
- shutil.move( os.path.join( html_path, 'fimo.xml' ), xml_out )
- shutil.move( os.path.join( html_path, 'fimo.html' ), html_out )
-
- out_file = open( interval_out, 'wb' )
- out_file.write( "#%s\n" % "\t".join( ( "chr", "start", "end", "pattern name", "score", "strand", "matched sequence", "p-value", "q-value" ) ) )
- for line in open( txt_out ):
- if line.startswith( '#' ): continue
- fields = line.rstrip( "\n\r" ).split( "\t" )
- start, end = int( fields[2] ), int( fields[3] )
- sequence = fields[7]
- if start > end:
- start, end = end, start #flip start and end, and set strand
- strand = "-"
- sequence = DNA_reverse_complement( sequence ) #we want sequences relative to strand; FIMO always provides + stranded sequence
- else:
- strand = "+"
- start -= 1 #make 0-based start position
- out_file.write( "%s\n" % "\t".join( [ fields[1], str( start ), str( end ), fields[0], fields[4], strand, sequence, fields[5], fields[6] ] ) )
- out_file.close()
- if __name__ == "__main__": main()