|
0
|
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
|
|
|
9 See below for options
|
|
|
10 """
|
|
|
11
|
|
|
12 import optparse, os, shutil, subprocess, sys, tempfile
|
|
|
13
|
|
|
14 def stop_err( msg ):
|
|
|
15 sys.stderr.write( '%s\n' % msg )
|
|
|
16 sys.exit()
|
|
|
17
|
|
|
18 def parseRefLoc( refLoc, refUID ):
|
|
|
19 for line in open( refLoc ):
|
|
|
20 if not line.startswith( '#' ):
|
|
|
21 fields = line.strip().split( '\t' )
|
|
|
22 if len( fields ) >= 3:
|
|
|
23 if fields[0] == refUID:
|
|
|
24 return fields[1]
|
|
|
25 return None
|
|
|
26
|
|
|
27 def __main__():
|
|
|
28 #Parse Command Line
|
|
|
29 parser = optparse.OptionParser()
|
|
|
30 parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to index and use' )
|
|
|
31 parser.add_option( '-u', '--refUID', dest='refUID', help='The pre-index reference genome unique Identifier' )
|
|
|
32 parser.add_option( '-L', '--refLocations', dest='refLocations', help='The filepath to the srma indices location file' )
|
|
|
33 parser.add_option( '-i', '--input', dest='input', help='The SAM/BAM input file' )
|
|
|
34 parser.add_option( '-I', '--inputIndex', dest='inputIndex', help='The SAM/BAM input index 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( '-j', '--jarBin', dest='jarBin', default='', help='The path to where jars are stored' )
|
|
|
47 parser.add_option( '-f', '--jarFile', dest='jarFile', help='The file name of the jar file to use')
|
|
|
48 (options, args) = parser.parse_args()
|
|
|
49
|
|
|
50 # make temp directory for srma
|
|
|
51 tmp_dir = tempfile.mkdtemp()
|
|
|
52 buffsize = 1048576
|
|
|
53
|
|
|
54 # set up reference filenames
|
|
|
55 reference_filepath_name = None
|
|
|
56 # need to create SRMA dict and Samtools fai files for custom genome
|
|
|
57 if options.fileSource == 'history':
|
|
|
58 try:
|
|
|
59 reference_filepath = tempfile.NamedTemporaryFile( dir=tmp_dir, suffix='.fa' )
|
|
|
60 reference_filepath_name = reference_filepath.name
|
|
|
61 reference_filepath.close()
|
|
|
62 fai_filepath_name = '%s.fai' % reference_filepath_name
|
|
|
63 dict_filepath_name = reference_filepath_name.replace( '.fa', '.dict' )
|
|
|
64 os.symlink( options.ref, reference_filepath_name )
|
|
|
65 # create fai file using Samtools
|
|
|
66 index_fai_cmd = 'samtools faidx %s' % reference_filepath_name
|
|
|
67 try:
|
|
|
68 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
69 tmp_stderr = open( tmp, 'wb' )
|
|
|
70 proc = subprocess.Popen( args=index_fai_cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
71 returncode = proc.wait()
|
|
|
72 tmp_stderr.close()
|
|
|
73 # get stderr, allowing for case where it's very large
|
|
|
74 tmp_stderr = open( tmp, 'rb' )
|
|
|
75 stderr = ''
|
|
|
76 try:
|
|
|
77 while True:
|
|
|
78 stderr += tmp_stderr.read( buffsize )
|
|
|
79 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
80 break
|
|
|
81 except OverflowError:
|
|
|
82 pass
|
|
|
83 tmp_stderr.close()
|
|
|
84 if returncode != 0:
|
|
|
85 raise Exception, stderr
|
|
|
86 except Exception, e:
|
|
|
87 # clean up temp dir
|
|
|
88 if os.path.exists( tmp_dir ):
|
|
|
89 shutil.rmtree( tmp_dir )
|
|
|
90 stop_err( 'Error creating Samtools index for custom genome file: %s\n' % str( e ) )
|
|
|
91 # create dict file using SRMA
|
|
|
92 dict_cmd = 'java -cp "%s" net.sf.picard.sam.CreateSequenceDictionary R=%s O=%s' % ( os.path.join( options.jarBin, options.jarFile ), reference_filepath_name, dict_filepath_name )
|
|
|
93 try:
|
|
|
94 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
95 tmp_stderr = open( tmp, 'wb' )
|
|
|
96 proc = subprocess.Popen( args=dict_cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
97 returncode = proc.wait()
|
|
|
98 tmp_stderr.close()
|
|
|
99 # get stderr, allowing for case where it's very large
|
|
|
100 tmp_stderr = open( tmp, 'rb' )
|
|
|
101 stderr = ''
|
|
|
102 try:
|
|
|
103 while True:
|
|
|
104 stderr += tmp_stderr.read( buffsize )
|
|
|
105 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
106 break
|
|
|
107 except OverflowError:
|
|
|
108 pass
|
|
|
109 tmp_stderr.close()
|
|
|
110 if returncode != 0:
|
|
|
111 raise Exception, stderr
|
|
|
112 except Exception, e:
|
|
|
113 # clean up temp dir
|
|
|
114 if os.path.exists( tmp_dir ):
|
|
|
115 shutil.rmtree( tmp_dir )
|
|
|
116 stop_err( 'Error creating index for custom genome file: %s\n' % str( e ) )
|
|
|
117 except Exception, e:
|
|
|
118 # clean up temp dir
|
|
|
119 if os.path.exists( tmp_dir ):
|
|
|
120 shutil.rmtree( tmp_dir )
|
|
|
121 stop_err( 'Problem handling SRMA index (dict file) for custom genome file: %s\n' % str( e ) )
|
|
|
122 # using built-in dict/index files
|
|
|
123 else:
|
|
|
124 if options.ref:
|
|
|
125 reference_filepath_name = options.ref
|
|
|
126 else:
|
|
|
127 reference_filepath_name = parseRefLoc( options.refLocation, options.refUID )
|
|
|
128 if reference_filepath_name is None:
|
|
|
129 raise ValueError( 'A valid genome reference was not provided.' )
|
|
|
130
|
|
|
131 # set up aligning and generate aligning command options
|
|
|
132 if options.params == 'pre_set':
|
|
|
133 srma_cmds = ''
|
|
|
134 else:
|
|
|
135 if options.useSequenceQualities == 'true':
|
|
|
136 useSequenceQualities = 'true'
|
|
|
137 else:
|
|
|
138 useSequenceQualities = 'false'
|
|
|
139 ranges = 'null'
|
|
|
140 if options.range == 'None':
|
|
|
141 range = 'null'
|
|
|
142 else:
|
|
|
143 range = options.range
|
|
|
144 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 )
|
|
|
145
|
|
|
146 # perform alignments
|
|
|
147 buffsize = 1048576
|
|
|
148 try:
|
|
|
149 #symlink input bam and index files due to the naming conventions required by srma here
|
|
|
150 input_bam_filename = os.path.join( tmp_dir, '%s.bam' % os.path.split( options.input )[-1] )
|
|
|
151 os.symlink( options.input, input_bam_filename )
|
|
|
152 input_bai_filename = "%s.bai" % os.path.splitext( input_bam_filename )[0]
|
|
|
153 os.symlink( options.inputIndex, input_bai_filename )
|
|
|
154
|
|
|
155 #create a temp output name, ending in .bam due to required naming conventions? unkown if required
|
|
|
156 output_bam_filename = os.path.join( tmp_dir, "%s.bam" % os.path.split( options.output )[-1] )
|
|
|
157 # generate commandline
|
|
|
158 cmd = 'java -jar %s I=%s O=%s R=%s %s' % ( os.path.join( options.jarBin, options.jarFile ), input_bam_filename, output_bam_filename, reference_filepath_name, srma_cmds )
|
|
|
159
|
|
|
160 # need to nest try-except in try-finally to handle 2.4
|
|
|
161 try:
|
|
|
162 try:
|
|
|
163 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
164 tmp_stderr = open( tmp, 'wb' )
|
|
|
165 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
166 returncode = proc.wait()
|
|
|
167 tmp_stderr.close()
|
|
|
168 # get stderr, allowing for case where it's very large
|
|
|
169 tmp_stderr = open( tmp, 'rb' )
|
|
|
170 stderr = ''
|
|
|
171 try:
|
|
|
172 while True:
|
|
|
173 stderr += tmp_stderr.read( buffsize )
|
|
|
174 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
175 break
|
|
|
176 except OverflowError:
|
|
|
177 pass
|
|
|
178 tmp_stderr.close()
|
|
|
179 if returncode != 0:
|
|
|
180 raise Exception, stderr
|
|
|
181 except Exception, e:
|
|
|
182 raise Exception, 'Error executing SRMA. ' + str( e )
|
|
|
183 # move file from temp location (with .bam name) to provided path
|
|
|
184 shutil.move( output_bam_filename, options.output )
|
|
|
185 # check that there are results in the output file
|
|
|
186 if os.path.getsize( options.output ) <= 0:
|
|
|
187 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.'
|
|
|
188 except Exception, e:
|
|
|
189 stop_err( 'The re-alignment failed.\n' + str( e ) )
|
|
|
190 finally:
|
|
|
191 # clean up temp dir
|
|
|
192 if os.path.exists( tmp_dir ):
|
|
|
193 shutil.rmtree( tmp_dir )
|
|
|
194
|
|
|
195 if __name__=="__main__": __main__()
|
