/tools/meme/fimo_wrapper.py

 1#!/usr/bin/env python
 2#Dan Blankenberg
 3
 4"""
 5Read text output from FIMO and create an interval file.
 6"""
 7import sys, tempfile, subprocess, shutil, os
 8from galaxy_utils.sequence.transform import DNA_reverse_complement
 9
10buffsize = 1048576
11
12def stop_err( msg ):
13    sys.stderr.write( msg )
14    sys.exit()
15
16def main():
17    assert len( sys.argv ) == 8, "Wrong number of arguments"
18    sys.argv.pop(0)
19    fimo_cmd = sys.argv.pop(0)
20    html_path = sys.argv.pop(0)
21    html_out = sys.argv.pop(0)
22    interval_out = sys.argv.pop(0)
23    txt_out = sys.argv.pop(0)
24    xml_out = sys.argv.pop(0)
25    gff_out = sys.argv.pop(0)
26    
27    #run fimo
28    try:
29        tmp_stderr = tempfile.NamedTemporaryFile()
30        #tmp_stderr = open( tmp_filename, 'wb' )
31        proc = subprocess.Popen( args=fimo_cmd, shell=True, stderr=tmp_stderr )
32        returncode = proc.wait()
33        #tmp_stderr.close()
34        # get stderr, allowing for case where it's very large
35        #tmp_stderr = open( tmp, 'rb' )
36        tmp_stderr.seek(0)
37        stderr = ''
38        try:
39            while True:
40                stderr += tmp_stderr.read( buffsize )
41                if not stderr or len( stderr ) % buffsize != 0:
42                    break
43        except OverflowError:
44            pass
45        
46        if returncode != 0:
47            raise Exception, stderr
48    except Exception, e:
49        raise Exception, 'Error running FIMO:\n' + str( e )
50
51    shutil.move( os.path.join( html_path, 'fimo.txt' ), txt_out )
52    shutil.move( os.path.join( html_path, 'fimo.gff' ), gff_out )
53    shutil.move( os.path.join( html_path, 'fimo.xml' ), xml_out )
54    shutil.move( os.path.join( html_path, 'fimo.html' ), html_out )
55    
56    out_file = open( interval_out, 'wb' )
57    out_file.write( "#%s\n" % "\t".join( ( "chr", "start", "end", "pattern name", "score", "strand", "matched sequence", "p-value", "q-value" ) ) )
58    for line in open( txt_out ):
59        if line.startswith( '#' ): continue
60        fields = line.rstrip( "\n\r" ).split( "\t" )
61        start, end = int( fields[2] ), int( fields[3] )
62        sequence = fields[7]
63        if start > end:
64            start, end = end, start #flip start and end, and set strand
65            strand = "-"
66            sequence = DNA_reverse_complement( sequence ) #we want sequences relative to strand; FIMO always provides + stranded sequence
67        else:
68            strand = "+"
69        start -= 1 #make 0-based start position
70        out_file.write( "%s\n" % "\t".join( [ fields[1], str( start ), str( end ), fields[0], fields[4], strand, sequence, fields[5], fields[6] ] ) )
71    out_file.close()
72
73if __name__ == "__main__": main()