annotate SNV/annotate.py @ 7:351b3acadd17 default tip

Uploaded

#!/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__()