Mercurial > repos > ryanmorin > nextgen_variant_identification
view SNV/annotate.py @ 7:351b3acadd17 default tip
Uploaded
author | ryanmorin |
---|---|
date | Tue, 18 Oct 2011 18:33:15 -0400 |
parents | 74f5ea818cea |
children |
line wrap: on
line source
#!/usr/bin/env python """ 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) usage: %prog [options] -i, --input1=i: bam file -o, --output1=o: Output pileup """ import os, shutil, subprocess, sys, tempfile from galaxy import eggs import pkg_resources; pkg_resources.require( "bx-python" ) from bx.cookbook import doc_optparse def stop_err( msg ): sys.stderr.write( '%s\n' % msg ) sys.exit() def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ): seqFile = '%s/sam_fa_indices.loc' % GALAXY_DATA_INDEX_DIR seqPath = '' for line in open( seqFile ): line = line.rstrip( '\r\n' ) if line and not line.startswith( '#' ) and line.startswith( 'index' ): fields = line.split( '\t' ) if len( fields ) < 3: continue if fields[1] == dbkey: seqPath = fields[2].strip() break return seqPath def __main__(): #Parse Command Line options, args = doc_optparse.parse( __doc__ ) #make temp dir tmpDir = tempfile.mkdtemp() cmd = 'identify_nonsynonymous_mutations.pl < %s > %s' try: # have to nest try-except in try-finally to handle 2.4 try: cmd = cmd % ( options.input1, options.output1 ) print(cmd, "\n") tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name tmp_stderr = open( tmp, 'wb' ) proc = subprocess.Popen( args=cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() ) returncode = proc.wait() tmp_stderr.close() #did it succeed? # get stderr, allowing for case where it's very large tmp_stderr = open( tmp, 'rb' ) stderr = '' try: while True: stderr += tmp_stderr.read( ) if not stderr: break except OverflowError: pass tmp_stderr.close() if returncode != 0: raise Exception, stderr except Exception, e: stop_err( 'Error running annotator tool\n' + str( e ) ) finally: #clean up temp files if os.path.exists( tmpDir ): shutil.rmtree( tmpDir ) # check that there are results in the output file if os.path.getsize( options.output1 ) > 0: sys.stdout.write( 'wrote annotated output' ) else: stop_err( 'The output file is empty' ) if __name__ == "__main__" : __main__()