Mercurial > repos > cjav > smalt

comparison tools/smalt_wrapper.py @ 0:747433a6de00 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit.
author cjav
date Wed, 13 Feb 2013 13:27:44 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:747433a6de00
1 #!/usr/bin/env python
2
3 """
4 Runs Smalt on single-end or paired-end data.
5 Produces a SAM file containing the mappings.
6 Works with Smalt version 0.7.1.
7
8 usage: smalt_wrapper.py [options]
9
10 See below for options
11 """
12
13 import optparse, os, shutil, subprocess, sys, tempfile
14
15 def stop_err( msg ):
16 sys.stderr.write( '%s\n' % msg )
17 sys.exit()
18
19 def __main__():
20 #Parse Command Line
21 parser = optparse.OptionParser()
22 parser.add_option( '-t', '--threads', dest='threads', help='The number of threads to use' )
23 parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to use or index' )
24 parser.add_option( '-f', '--input1', dest='fastq', help='The (forward) fastq file to use for the mapping' )
25 parser.add_option( '-F', '--input2', dest='rfastq', help='The reverse fastq file to use for mapping if paired-end data' )
26 parser.add_option( '-u', '--output', dest='output', help='The file to save the output (SAM format)' )
27 parser.add_option( '-g', '--genAlignType', dest='genAlignType', help='The type of pairing (single or paired)' )
28 parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
29 parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one form history (indexed or history)' )
30 parser.add_option( '-x', '--exhaustiveSearch', dest='exhaustiveSearch', help='This flag triggers a more exhaustive search for alignments at the cost of decreased speed' )
31 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' )
32 parser.add_option( '-d', '--scorDiff', dest='scorDiff', help='Set a threshold of the Smith-Waterman alignment score relative to the maximum score' )
33 parser.add_option( '-i', '--insertMax', dest='insertMax', help='Maximum insert size (Only in paired-end mode)' )
34 parser.add_option( '-j', '--insertMin', dest='insertMin', help='Minimum insert size (Only in paired-end mode)' )
35 parser.add_option( '-l', '--pairTyp', dest='pairTyp', help='Type of read pair library, can be either pe, mp or pp' )
36 parser.add_option( '-m', '--minScor', dest='minScor', help='Sets an absolute threshold of the Smith-Waterman scores' )
37 parser.add_option( '-a', '--partialAlignments', dest='partialAlignments', help='Report partial alignments if they are complementary on the read (split reads)' )
38 parser.add_option( '-q', '--minBasq', dest='minBasq', help='Sets a base quality threshold (0 <= minbasq <= 10, default 0)' )
39 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".' )
40 parser.add_option( '-w', '--complexityWeighted', dest='complexityWeighted', help='Smith-Waterman scores are complexity weighted' )
41 parser.add_option( '-y', '--minId', dest='minId', help='Sets an identity threshold for a mapping to be reported' )
42 parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
43 parser.add_option( '-X', '--do_not_build_index', dest='do_not_build_index', action='store_true', help="Don't build index" )
44 parser.add_option( '-H', '--suppressHeader', dest='suppressHeader', help='Suppress header' )
45 (options, args) = parser.parse_args()
46
47 # output version # of tool
48 try:
49 tmp = tempfile.NamedTemporaryFile().name
50 tmp_stdout = open( tmp, 'wb' )
51 proc = subprocess.Popen( args='smalt 2>&1', shell=True, stdout=tmp_stdout )
52 tmp_stdout.close()
53 returncode = proc.wait()
54 stdout = None
55 for line in open( tmp_stdout.name, 'rb' ):
56 if line.lower().find( 'version' ) >= 0:
57 stdout = line.strip()
58 break
59 if stdout:
60 sys.stdout.write( 'SMALT %s\n' % stdout )
61 else:
62 raise Exception
63 except:
64 sys.stdout.write( 'Could not determine SMALT version\n' )
65
66 fastq = options.fastq
67 if options.rfastq:
68 rfastq = options.rfastq
69
70 # make temp directory for placement of indices
71 tmp_index_dir = tempfile.mkdtemp()
72 tmp_dir = tempfile.mkdtemp()
73 # index if necessary
74 if options.fileSource == 'history' and not options.do_not_build_index:
75 ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
76 ref_file_name = ref_file.name
77 ref_file.close()
78 os.symlink( options.ref, ref_file_name )
79 cmd1 = 'smalt index %s %s' % ( ref_file_name, ref_file_name )
80 try:
81 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
82 tmp_stderr = open( tmp, 'wb' )
83 proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_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 buffsize = 1048576
90 try:
91 while True:
92 stderr += tmp_stderr.read( buffsize )
93 if not stderr or len( stderr ) % buffsize != 0:
94 break
95 except OverflowError:
96 pass
97 tmp_stderr.close()
98 if returncode != 0:
99 raise Exception, stderr
100 except Exception, e:
101 # clean up temp dirs
102 if os.path.exists( tmp_index_dir ):
103 shutil.rmtree( tmp_index_dir )
104 if os.path.exists( tmp_dir ):
105 shutil.rmtree( tmp_dir )
106 stop_err( 'Error indexing reference sequence. ' + str( e ) )
107 else:
108 ref_file_name = options.ref
109
110 # set up aligning and generate aligning command options
111 if options.params == 'pre_set':
112 aligning_cmds = '-n %s ' % ( options.threads )
113 gen_alignment_cmds = ''
114 else:
115 if options.exhaustiveSearch == 'true':
116 exhaustiveSearch = '-x'
117 minCover = '-c %s' % options.minCover
118 else:
119 exhaustiveSearch = ''
120 minCover = ''
121 if options.partialAlignments == 'true':
122 partialAlignments = '-x'
123 else:
124 partialAlignments = ''
125 if options.complexityWeighted == 'true':
126 complexityWeighted = '-w'
127 else:
128 complexityWeighted = ''
129 aligning_cmds = '-d %s -m %s -q %s -r %s -y %s' % \
130 ( options.scorDiff, options.minScor, options.minBasq, options.seed, options.minId )
131 if options.genAlignType == 'paired':
132 gen_alignment_cmds = '-i %s -j %s -l %s' % ( options.insertMax, options.insertMin, options.pairTyp )
133 else:
134 gen_alignment_cmds = ''
135 # prepare actual aligning and generate aligning commands
136 if options.genAlignType == 'paired':
137 cmd = 'smalt map %s %s -o %s %s %s ' % ( aligning_cmds, gen_alignment_cmds, options.output, ref_file_name, fastq, rfastq )
138 else:
139 cmd = 'smalt map %s -o %s %s %s ' % ( aligning_cmds, options.output, ref_file_name, fastq )
140 # perform alignments
141 buffsize = 1048576
142 try:
143 # need to nest try-except in try-finally to handle 2.4
144 try:
145 # align
146 try:
147 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
148 tmp_stderr = open( tmp, 'wb' )
149 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
150 returncode = proc.wait()
151 tmp_stderr.close()
152 # get stderr, allowing for case where it's very large
153 tmp_stderr = open( tmp, 'rb' )
154 stderr = ''
155 try:
156 while True:
157 stderr += tmp_stderr.read( buffsize )
158 if not stderr or len( stderr ) % buffsize != 0:
159 break
160 except OverflowError:
161 pass
162 tmp_stderr.close()
163 if returncode != 0:
164 raise Exception, stderr
165 except Exception, e:
166 raise Exception, 'Error aligning sequence. ' + str( e )
167 # remove header if necessary
168 if options.suppressHeader == 'true':
169 tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
170 tmp_out_name = tmp_out.name
171 tmp_out.close()
172 try:
173 shutil.move( options.output, tmp_out_name )
174 except Exception, e:
175 raise Exception, 'Error moving output file before removing headers. ' + str( e )
176 fout = file( options.output, 'w' )
177 for line in file( tmp_out.name, 'r' ):
178 if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
179 fout.write( line )
180 fout.close()
181 # check that there are results in the output file
182 if os.path.getsize( options.output ) > 0:
183 sys.stdout.write( 'SMALT run on %s-end data' % options.genAlignType )
184 else:
185 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.'
186 except Exception, e:
187 stop_err( 'The alignment failed.\n' + str( e ) )
188 finally:
189 # clean up temp dir
190 if os.path.exists( tmp_index_dir ):
191 shutil.rmtree( tmp_index_dir )
192 if os.path.exists( tmp_dir ):
193 shutil.rmtree( tmp_dir )
194
195 if __name__=="__main__": __main__()