Mercurial > repos > nilshomer > srma_wrapper

comparison srma_wrapper.py @ 0:9d60d2fce247 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
author nilshomer
date Tue, 07 Jun 2011 17:43:07 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9d60d2fce247
1 #!/usr/bin/env python
2
3 """
4 Runs SRMA on a SAM/BAM file;
5 TODO: more documentation
6
7 usage: srma_wrapper.py [options]
8 -r, --ref=r: The reference genome to use or index
9 -i, --input=i: The SAM/BAM input file
10 -o, --output=o: The SAM/BAM output file
11 -O, --offset=O: The alignment offset
12 -Q, --minMappingQuality=Q: The minimum mapping quality
13 -P, --minAlleleProbability=P: The minimum allele probability conditioned on coverage (for the binomial quantile).
14 -C, --minAlleleCoverage=C: The minimum haploid coverage for the consensus. Default value: 3. This option can be set
15 -R, --range=R: A range to examine
16 -c, --correctBases=c: Correct bases
17 -q, --useSequenceQualities=q: Use sequence qualities
18 -M, --maxHeapSize=M: The maximum number of nodes on the heap before re-alignment is ignored
19 -s, --fileSource=s: Whether to use a previously indexed reference sequence or one from history (indexed or history)
20 -p, --params=p: Parameter setting to use (pre_set or full)
21 -D, --dbkey=D: Dbkey for reference genome
22 """
23
24 import optparse, os, shutil, subprocess, sys, tempfile
25
26 def stop_err( msg ):
27 sys.stderr.write( '%s\n' % msg )
28 sys.exit()
29
30 def __main__():
31 #Parse Command Line
32 parser = optparse.OptionParser()
33 parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to use or index' )
34 parser.add_option( '-i', '--input', dest='input', help='The SAM/BAM input file' )
35 parser.add_option( '-o', '--output', dest='output', help='The SAM/BAM output file' )
36 parser.add_option( '-O', '--offset', dest='offset', help='The alignment offset' )
37 parser.add_option( '-Q', '--minMappingQuality', dest='minMappingQuality', help='The minimum mapping quality' )
38 parser.add_option( '-P', '--minAlleleProbability', dest='minAlleleProbability', help='The minimum allele probability conditioned on coverage (for the binomial quantile).' )
39 parser.add_option( '-C', '--minAlleleCoverage', dest='minAlleleCoverage', help='The minimum haploid coverage for the consensus' )
40 parser.add_option( '-R', '--range', dest='range', help='A range to examine' )
41 parser.add_option( '-c', '--correctBases', dest='correctBases', help='Correct bases ' )
42 parser.add_option( '-q', '--useSequenceQualities', dest='useSequenceQualities', help='Use sequence qualities ' )
43 parser.add_option( '-M', '--maxHeapSize', dest='maxHeapSize', help='The maximum number of nodes on the heap before re-alignment is ignored' )
44 parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one from history (indexed or history)' )
45 parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
46 parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
47 (options, args) = parser.parse_args()
48
49 # make temp directory for bfast
50 tmp_dir = '%s/' % tempfile.mkdtemp()
51
52 # assume indexing has already been done
53 if options.fileSource == 'history':
54 stop_err( 'Error: indexing not implemented' )
55
56 # set up aligning and generate aligning command options
57 if options.params == 'pre_set':
58 srma_cmds = ''
59 else:
60 if options.useSequenceQualities == 'true':
61 useSequenceQualities = 'true'
62 else:
63 useSequenceQualities = 'false'
64 ranges = 'null'
65 if options.range == 'None':
66 range = 'null'
67 else:
68 range = options.range
69
70 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 )
71
72 # how do we call a JAR file?
73 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 )
74
75 # perform alignments
76 buffsize = 1048576
77 try:
78 # need to nest try-except in try-finally to handle 2.4
79 try:
80 try:
81 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
82 tmp_stderr = open( tmp, 'wb' )
83 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
84 returncode = proc.wait()
85 tmp_stderr.close()
86 # get stderr, allowing for case where it's very large
87 tmp_stderr = open( tmp, 'rb' )
88 stderr = ''
89 try:
90 while True:
91 stderr += tmp_stderr.read( buffsize )
92 if not stderr or len( stderr ) % buffsize != 0:
93 break
94 except OverflowError:
95 pass
96 tmp_stderr.close()
97 if returncode != 0:
98 raise Exception, stderr
99 except Exception, e:
100 raise Exception, 'Error executing SRMA. ' + str( e )
101 # check that there are results in the output file
102 if os.path.getsize( options.output ) > 0:
103 if "0" == options.space:
104 sys.stdout.write( 'BFAST run on Base Space data' )
105 else:
106 sys.stdout.write( 'BFAST run on Color Space data' )
107 else:
108 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.'
109 except Exception, e:
110 stop_err( 'The alignment failed.\n' + str( e ) )
111 finally:
112 # clean up temp dir
113 if os.path.exists( tmp_dir ):
114 shutil.rmtree( tmp_dir )
115
116 if __name__=="__main__": __main__()