diff tools/meme/fimo_wrapper.py @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/meme/fimo_wrapper.py	Fri Mar 09 19:37:19 2012 -0500
@@ -0,0 +1,73 @@
+#!/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()