Mercurial > repos > nilshomer > srma_wrapper

annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9d60d2fce247 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
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
2
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
3 """
9d60d2fce247 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;
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
5 TODO: more documentation
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
6
9d60d2fce247 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]
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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).
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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
9d60d2fce247 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)
9d60d2fce247 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)
9d60d2fce247 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
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
22 """
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
23
9d60d2fce247 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
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
25
9d60d2fce247 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 ):
9d60d2fce247 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 )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
28 sys.exit()
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
29
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
30 def __main__():
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
31 #Parse Command Line
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
32 parser = optparse.OptionParser()
9d60d2fce247 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' )
9d60d2fce247 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' )
9d60d2fce247 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' )
9d60d2fce247 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' )
9d60d2fce247 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' )
9d60d2fce247 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).' )
9d60d2fce247 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' )
9d60d2fce247 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' )
9d60d2fce247 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 ' )
9d60d2fce247 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 ' )
9d60d2fce247 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' )
9d60d2fce247 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)' )
9d60d2fce247 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)' )
9d60d2fce247 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' )
9d60d2fce247 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()
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
48
9d60d2fce247 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
9d60d2fce247 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()
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
51
9d60d2fce247 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
9d60d2fce247 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':
9d60d2fce247 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' )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
55
9d60d2fce247 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
9d60d2fce247 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':
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
58 srma_cmds = ''
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
59 else:
9d60d2fce247 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':
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
61 useSequenceQualities = 'true'
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
62 else:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
63 useSequenceQualities = 'false'
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
64 ranges = 'null'
9d60d2fce247 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':
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
66 range = 'null'
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
67 else:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
68 range = options.range
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
69
9d60d2fce247 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 )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
71
9d60d2fce247 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?
9d60d2fce247 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 )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
74
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
75 # perform alignments
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
76 buffsize = 1048576
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
77 try:
9d60d2fce247 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
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
79 try:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
80 try:
9d60d2fce247 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
9d60d2fce247 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' )
9d60d2fce247 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() )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
84 returncode = proc.wait()
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
85 tmp_stderr.close()
9d60d2fce247 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
9d60d2fce247 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' )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
88 stderr = ''
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
89 try:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
90 while True:
9d60d2fce247 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 )
9d60d2fce247 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:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
93 break
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
94 except OverflowError:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
95 pass
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
96 tmp_stderr.close()
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
97 if returncode != 0:
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
98 raise Exception, stderr
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
99 except Exception, e:
9d60d2fce247 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 )
9d60d2fce247 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
9d60d2fce247 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:
9d60d2fce247 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:
9d60d2fce247 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' )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
105 else:
9d60d2fce247 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' )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
107 else:
9d60d2fce247 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.'
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
109 except Exception, e:
9d60d2fce247 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 ) )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
111 finally:
9d60d2fce247 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
9d60d2fce247 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 ):
9d60d2fce247 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 )
9d60d2fce247 Migrated tool version 0.1.1 from old tool shed archive to new tool shed repository
nilshomer
parents:
diff changeset
115
9d60d2fce247 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__()