Mercurial > repos > ryanmorin > nextgen_variant_identification

comparison SNV/annotate.py @ 0:74f5ea818cea

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author ryanmorin
date Wed, 12 Oct 2011 19:50:38 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:74f5ea818cea
1 #!/usr/bin/env python
2
3 """
4 Annotates filtered SNVMix output that has been attached to codon information (outputs the same data with additional columns describing the predicted effect of the SNV)
5
6 usage: %prog [options]
7 -i, --input1=i: bam file
8 -o, --output1=o: Output pileup
9 """
10
11 import os, shutil, subprocess, sys, tempfile
12 from galaxy import eggs
13 import pkg_resources; pkg_resources.require( "bx-python" )
14 from bx.cookbook import doc_optparse
15
16
17
18 def stop_err( msg ):
19 sys.stderr.write( '%s\n' % msg )
20 sys.exit()
21
22 def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ):
23 seqFile = '%s/sam_fa_indices.loc' % GALAXY_DATA_INDEX_DIR
24 seqPath = ''
25 for line in open( seqFile ):
26 line = line.rstrip( '\r\n' )
27 if line and not line.startswith( '#' ) and line.startswith( 'index' ):
28 fields = line.split( '\t' )
29 if len( fields ) < 3:
30 continue
31 if fields[1] == dbkey:
32 seqPath = fields[2].strip()
33 break
34 return seqPath
35
36 def __main__():
37 #Parse Command Line
38 options, args = doc_optparse.parse( __doc__ )
39 #make temp dir
40 tmpDir = tempfile.mkdtemp()
41
42 cmd = 'identify_nonsynonymous_mutations.pl < %s > %s'
43 try:
44 # have to nest try-except in try-finally to handle 2.4
45 try:
46 cmd = cmd % ( options.input1, options.output1 )
47 print(cmd, "\n")
48 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name
49 tmp_stderr = open( tmp, 'wb' )
50
51 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() )
52 returncode = proc.wait()
53 tmp_stderr.close()
54 #did it succeed?
55 # get stderr, allowing for case where it's very large
56 tmp_stderr = open( tmp, 'rb' )
57 stderr = ''
58 try:
59 while True:
60 stderr += tmp_stderr.read( )
61 if not stderr:
62 break
63 except OverflowError:
64 pass
65 tmp_stderr.close()
66 if returncode != 0:
67 raise Exception, stderr
68 except Exception, e:
69 stop_err( 'Error running annotator tool\n' + str( e ) )
70 finally:
71 #clean up temp files
72 if os.path.exists( tmpDir ):
73 shutil.rmtree( tmpDir )
74 # check that there are results in the output file
75 if os.path.getsize( options.output1 ) > 0:
76 sys.stdout.write( 'wrote annotated output' )
77 else:
78 stop_err( 'The output file is empty' )
79
80 if __name__ == "__main__" : __main__()