diff tools/smalt_wrapper.py @ 1:54855bd8d107 draft default tip

First attempt to get tool_dependencies.xml right.
author cjav
date Wed, 20 Mar 2013 17:07:14 -0400
parents 747433a6de00
children
line wrap: on
line diff
--- a/tools/smalt_wrapper.py	Wed Feb 13 13:27:44 2013 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,195 +0,0 @@
-#!/usr/bin/env python
-
-"""
-Runs Smalt on single-end or paired-end data.
-Produces a SAM file containing the mappings.
-Works with Smalt version 0.7.1.
-
-usage: smalt_wrapper.py [options]
-
-See below for options
-"""
-
-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( '-t', '--threads', dest='threads', help='The number of threads to use' )
-    parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to use or index' )
-    parser.add_option( '-f', '--input1', dest='fastq', help='The (forward) fastq file to use for the mapping' )
-    parser.add_option( '-F', '--input2', dest='rfastq', help='The reverse fastq file to use for mapping if paired-end data' )
-    parser.add_option( '-u', '--output', dest='output', help='The file to save the output (SAM format)' )
-    parser.add_option( '-g', '--genAlignType', dest='genAlignType', help='The type of pairing (single or paired)' )
-    parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
-    parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one form history (indexed or history)' )
-    parser.add_option( '-x', '--exhaustiveSearch', dest='exhaustiveSearch', help='This flag triggers a more exhaustive search for alignments at the cost of decreased speed' )
-    parser.add_option( '-c', '--minCover', dest='minCover', help='Only consider mappings where the k-mer word seeds cover the query read to a minimum extent' )
-    parser.add_option( '-d', '--scorDiff', dest='scorDiff', help='Set a threshold of the Smith-Waterman alignment score relative to the maximum score' )
-    parser.add_option( '-i', '--insertMax', dest='insertMax', help='Maximum insert size (Only in paired-end mode)' )
-    parser.add_option( '-j', '--insertMin', dest='insertMin', help='Minimum insert size (Only in paired-end mode)' )
-    parser.add_option( '-l', '--pairTyp', dest='pairTyp', help='Type of read pair library, can be either pe, mp or pp' )
-    parser.add_option( '-m', '--minScor', dest='minScor', help='Sets an absolute threshold of the Smith-Waterman scores' )
-    parser.add_option( '-a', '--partialAlignments', dest='partialAlignments', help='Report partial alignments if they are complementary on the read (split reads)' )
-    parser.add_option( '-q', '--minBasq', dest='minBasq', help='Sets a base quality threshold (0 <= minbasq <= 10, default 0)' )
-    parser.add_option( '-e', '--seed', dest='seed', help='If <seed> >= 0 report an alignment selected at random where there are multiple mappings with the same best alignment score. With <seed> = 0 (default) a seed is derived from the current calendar time. If <seed> < 0 reads with multiple best mappings are reported as "not mapped".' )
-    parser.add_option( '-w', '--complexityWeighted', dest='complexityWeighted', help='Smith-Waterman scores are complexity weighted' )
-    parser.add_option( '-y', '--minId', dest='minId', help='Sets an identity threshold for a mapping to be reported' )
-    parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
-    parser.add_option( '-X', '--do_not_build_index', dest='do_not_build_index', action='store_true', help="Don't build index" )
-    parser.add_option( '-H', '--suppressHeader', dest='suppressHeader', help='Suppress header' )
-    (options, args) = parser.parse_args()
-
-    # output version # of tool
-    try:
-        tmp = tempfile.NamedTemporaryFile().name
-        tmp_stdout = open( tmp, 'wb' )
-        proc = subprocess.Popen( args='smalt 2>&1', shell=True, stdout=tmp_stdout )
-        tmp_stdout.close()
-        returncode = proc.wait()
-        stdout = None
-        for line in open( tmp_stdout.name, 'rb' ):
-            if line.lower().find( 'version' ) >= 0:
-                stdout = line.strip()
-                break
-        if stdout:
-            sys.stdout.write( 'SMALT %s\n' % stdout )
-        else:
-            raise Exception
-    except:
-        sys.stdout.write( 'Could not determine SMALT version\n' )
-
-    fastq = options.fastq
-    if options.rfastq:
-         rfastq = options.rfastq
-
-    # make temp directory for placement of indices
-    tmp_index_dir = tempfile.mkdtemp()
-    tmp_dir = tempfile.mkdtemp()
-    # index if necessary
-    if options.fileSource == 'history' and not options.do_not_build_index:
-        ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
-        ref_file_name = ref_file.name
-        ref_file.close()
-        os.symlink( options.ref, ref_file_name )
-        cmd1 = 'smalt index %s %s' % ( ref_file_name, ref_file_name )
-        try:
-            tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
-            tmp_stderr = open( tmp, 'wb' )
-            proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_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 = ''
-            buffsize = 1048576
-            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:
-            # clean up temp dirs
-            if os.path.exists( tmp_index_dir ):
-                shutil.rmtree( tmp_index_dir )
-            if os.path.exists( tmp_dir ):
-                shutil.rmtree( tmp_dir )
-            stop_err( 'Error indexing reference sequence. ' + str( e ) )
-    else:
-        ref_file_name = options.ref
-
-    # set up aligning and generate aligning command options
-    if options.params == 'pre_set':
-        aligning_cmds = '-n %s ' % ( options.threads )
-        gen_alignment_cmds = ''
-    else:
-        if options.exhaustiveSearch == 'true':
-            exhaustiveSearch = '-x'
-            minCover = '-c %s' % options.minCover
-        else:
-            exhaustiveSearch = ''
-            minCover = ''
-        if options.partialAlignments == 'true':
-            partialAlignments = '-x'
-        else:
-            partialAlignments = ''
-        if options.complexityWeighted == 'true':
-            complexityWeighted = '-w'
-        else:
-            complexityWeighted = ''
-        aligning_cmds = '-d %s -m %s -q %s -r %s -y %s' % \
-                        ( options.scorDiff, options.minScor, options.minBasq, options.seed, options.minId )
-        if options.genAlignType == 'paired':
-            gen_alignment_cmds = '-i %s -j %s -l %s' % ( options.insertMax, options.insertMin, options.pairTyp )
-        else:
-            gen_alignment_cmds = ''
-    # prepare actual aligning and generate aligning commands
-    if options.genAlignType == 'paired':
-        cmd = 'smalt map %s %s -o %s %s %s ' % ( aligning_cmds, gen_alignment_cmds, options.output, ref_file_name, fastq, rfastq )
-    else:
-        cmd = 'smalt map %s -o %s %s %s ' % ( aligning_cmds, options.output, ref_file_name, fastq )
-    # perform alignments
-    buffsize = 1048576
-    try:
-        # need to nest try-except in try-finally to handle 2.4
-        try:
-            # align
-            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 aligning sequence. ' + str( e )
-            # remove header if necessary
-            if options.suppressHeader == 'true':
-                tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
-                tmp_out_name = tmp_out.name
-                tmp_out.close()
-                try:
-                    shutil.move( options.output, tmp_out_name )
-                except Exception, e:
-                    raise Exception, 'Error moving output file before removing headers. ' + str( e )
-                fout = file( options.output, 'w' )
-                for line in file( tmp_out.name, 'r' ):
-                    if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
-                        fout.write( line )
-                fout.close()
-            # check that there are results in the output file
-            if os.path.getsize( options.output ) > 0:
-                sys.stdout.write( 'SMALT run on %s-end data' % options.genAlignType )
-            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_index_dir ):
-            shutil.rmtree( tmp_index_dir )
-        if os.path.exists( tmp_dir ):
-            shutil.rmtree( tmp_dir )
-
-if __name__=="__main__": __main__()