Mercurial > repos > nilshomer > srma_wrapper

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