view srma_wrapper.py @ 0:9d60d2fce247 default tip

Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository

#!/usr/bin/env python

"""
Runs SRMA on a SAM/BAM file;
TODO: more documentation

usage: srma_wrapper.py [options]
    -r, --ref=r: The reference genome to use or index
    -i, --input=i: The SAM/BAM input file
    -o, --output=o: The SAM/BAM output file
    -O, --offset=O: The alignment offset
    -Q, --minMappingQuality=Q: The minimum mapping quality
    -P, --minAlleleProbability=P: The minimum allele probability conditioned on coverage (for the binomial quantile).
    -C, --minAlleleCoverage=C: The minimum haploid coverage for the consensus. Default value: 3. This option can be set 
    -R, --range=R: A range to examine
    -c, --correctBases=c: Correct bases 
    -q, --useSequenceQualities=q: Use sequence qualities 
    -M, --maxHeapSize=M: The maximum number of nodes on the heap before re-alignment is ignored
    -s, --fileSource=s: Whether to use a previously indexed reference sequence or one from history (indexed or history)
    -p, --params=p: Parameter setting to use (pre_set or full)
    -D, --dbkey=D: Dbkey for reference genome
"""

import optparse, os, shutil, subprocess, sys, tempfile

def stop_err( msg ):
    sys.stderr.write( '%s\n' % msg )
    sys.exit()

def __main__():
    #Parse Command Line
    parser = optparse.OptionParser()
    parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to use or index' )
    parser.add_option( '-i', '--input', dest='input', help='The SAM/BAM input file' )
    parser.add_option( '-o', '--output', dest='output', help='The SAM/BAM output file' )
    parser.add_option( '-O', '--offset', dest='offset', help='The alignment offset' )
    parser.add_option( '-Q', '--minMappingQuality', dest='minMappingQuality', help='The minimum mapping quality' )
    parser.add_option( '-P', '--minAlleleProbability', dest='minAlleleProbability', help='The minimum allele probability conditioned on coverage (for the binomial quantile).' )
    parser.add_option( '-C', '--minAlleleCoverage', dest='minAlleleCoverage', help='The minimum haploid coverage for the consensus' )
    parser.add_option( '-R', '--range', dest='range', help='A range to examine' )
    parser.add_option( '-c', '--correctBases', dest='correctBases', help='Correct bases ' )
    parser.add_option( '-q', '--useSequenceQualities', dest='useSequenceQualities', help='Use sequence qualities ' )
    parser.add_option( '-M', '--maxHeapSize', dest='maxHeapSize', help='The maximum number of nodes on the heap before re-alignment is ignored' )
    parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one from history (indexed or history)' )
    parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
    parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
    (options, args) = parser.parse_args()

    # make temp directory for bfast
    tmp_dir = '%s/' % tempfile.mkdtemp()
        
    # assume indexing has already been done
    if options.fileSource == 'history':
        stop_err( 'Error: indexing not implemented' )

    # set up aligning and generate aligning command options
    if options.params == 'pre_set':
        srma_cmds = ''
    else:
        if options.useSequenceQualities == 'true':
            useSequenceQualities = 'true'
        else:
            useSequenceQualities = 'false'
        ranges = 'null'
        if options.range == 'None':
            range = 'null'
        else:
            range = options.range
        
        srma_cmds = "OFFSET=%s MIN_MAPQ=%s MINIMUM_ALLELE_PROBABILITY=%s MINIMUM_ALLELE_COVERAGE=%s RANGES=%s RANGE=%s CORRECT_BASES=%s USE_SEQUENCE_QUALITIES=%s MAX_HEAP_SIZE=%s" % ( options.offset, options.minMappingQuality, options.minAlleleProbability, options.minAlleleCoverage, ranges, range, options.correctBases, options.useSequenceQualities, options.maxHeapSize )

    # how do we call a JAR file?
    cmd = 'java -jar /Users/nhomer/srma/build/jar/srma.jar I=%s O=%s R=%s %s' % ( options.input, options.output, options.ref, srma_cmds )
    
    # perform alignments
    buffsize = 1048576
    try:
        # need to nest try-except in try-finally to handle 2.4
        try:
            try:
                tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
                tmp_stderr = open( tmp, 'wb' )
                proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
                returncode = proc.wait()
                tmp_stderr.close()
                # get stderr, allowing for case where it's very large
                tmp_stderr = open( tmp, 'rb' )
                stderr = ''
                try:
                    while True:
                        stderr += tmp_stderr.read( buffsize )
                        if not stderr or len( stderr ) % buffsize != 0:
                            break
                except OverflowError:
                    pass
                tmp_stderr.close()
                if returncode != 0:
                    raise Exception, stderr
            except Exception, e:
                raise Exception, 'Error executing SRMA. ' + str( e )
            # check that there are results in the output file
            if os.path.getsize( options.output ) > 0:
                if "0" == options.space:
                    sys.stdout.write( 'BFAST run on Base Space data' )
                else:
                    sys.stdout.write( 'BFAST run on Color Space data' )
            else:
                raise Exception, 'The output file is empty. You may simply have no matches, or there may be an error with your input file or settings.'
        except Exception, e:
            stop_err( 'The alignment failed.\n' + str( e ) )
    finally:
        # clean up temp dir
        if os.path.exists( tmp_dir ):
            shutil.rmtree( tmp_dir )

if __name__=="__main__": __main__()