annotate bismark_wrapper.py @ 4:427fb56f2e41 draft default tip

- new options - fixes
author bjoern-gruening
date Fri, 01 Mar 2013 13:39:22 -0500
parents 36d124f44c0a
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
1 #!/usr/bin/env python
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
2
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
3 import argparse, os, shutil, subprocess, sys, tempfile, fileinput
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
4 import fileinput
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
5 from glob import glob
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
6
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
7 def stop_err( msg ):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
8 sys.stderr.write( "%s\n" % msg )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
9 sys.exit()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
10
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
11 def __main__():
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
12 #Parse Command Line
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
13 parser = argparse.ArgumentParser(description='Wrapper for the bismark bisulfite mapper.')
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
14 parser.add_argument( '-p', '--num-threads', dest='num_threads',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
15 type=int, default=4, help='Use this many threads to align reads. The default is 4.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
16
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
17 parser.add_argument( '--bismark_path', dest='bismark_path', help='Path to the bismark perl scripts' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
18
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
19 parser.add_argument( '--bowtie2', action='store_true', default=False, help='Running bismark with bowtie2 and not with bowtie.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
20
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
21 # input options
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
22 parser.add_argument( '--own-file', dest='own_file', help='' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
23 parser.add_argument( '-D', '--indexes-path', dest='index_path', help='Indexes directory; location of .ebwt and .fa files.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
24 parser.add_argument( '-O', '--output', dest='output' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
25 parser.add_argument( '--output-report-file', dest='output_report_file' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
26 parser.add_argument( '--suppress-header', dest='suppress_header', action="store_true" )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
27
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
28 parser.add_argument( '--mate-paired', dest='mate_paired', action='store_true', help='Reads are mate-paired', default=False)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
29
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
30
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
31 parser.add_argument( '-1', '--mate1', dest='mate1',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
32 help='The forward reads file in Sanger FASTQ or FASTA format.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
33 parser.add_argument( '-2', '--mate2', dest='mate2',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
34 help='The reverse reads file in Sanger FASTQ or FASTA format.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
35
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
36 parser.add_argument( '--output-unmapped-reads', dest='output_unmapped_reads',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
37 help='Additional output file with unmapped reads (single-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
38 parser.add_argument( '--output-unmapped-reads-l', dest='output_unmapped_reads_l',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
39 help='File name for unmapped reads (left, paired-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
40 parser.add_argument( '--output-unmapped-reads-r', dest='output_unmapped_reads_r',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
41 help='File name for unmapped reads (right, paired-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
42
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
43
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
44 parser.add_argument( '--output-suppressed-reads', dest='output_suppressed_reads',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
45 help='Additional output file with suppressed reads (single-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
46 parser.add_argument( '--output-suppressed-reads-l', dest='output_suppressed_reads_l',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
47 help='File name for suppressed reads (left, paired-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
48 parser.add_argument( '--output-suppressed-reads-r', dest='output_suppressed_reads_r',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
49 help='File name for suppressed reads (right, paired-end).' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
50
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
51
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
52 parser.add_argument( '--single-paired', dest='single_paired',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
53 help='The single-end reads file in Sanger FASTQ or FASTA format.' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
54
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
55 parser.add_argument( '--fastq', action='store_true', help='Query filetype is in FASTQ format')
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
56 parser.add_argument( '--fasta', action='store_true', help='Query filetype is in FASTA format')
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
57 parser.add_argument( '--phred64-quals', dest='phred64', action="store_true" )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
58
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
59
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
60 parser.add_argument( '--skip-reads', dest='skip_reads', type=int )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
61 parser.add_argument( '--qupto', type=int)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
62
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
63
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
64 # paired end options
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
65 parser.add_argument( '-I', '--minins', dest='min_insert' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
66 parser.add_argument( '-X', '--maxins', dest='max_insert' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
67 parser.add_argument( '--no-mixed', dest='no_mixed', action="store_true" )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
68 parser.add_argument( '--no-discordant', dest='no_discordant', action="store_true" )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
69
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
70 #parse general options
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
71 # default 20
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
72 parser.add_argument( '--seed-len', dest='seed_len', type=int)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
73 # default 15
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
74 parser.add_argument( '--seed-extention-attempts', dest='seed_extention_attempts', type=int )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
75 # default 0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
76 parser.add_argument( '--seed-mismatches', dest='seed_mismatches', type=int )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
77 # default 2
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
78 parser.add_argument( '--max-reseed', dest='max_reseed', type=int )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
79 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
80 # default 70
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
81 parser.add_argument( '--maqerr', dest='maqerr', type=int )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
82 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
83
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
84 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
85 The number of megabytes of memory a given thread is given to store path
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
86 descriptors in --best mode. Best-first search must keep track of many paths
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
87 at once to ensure it is always extending the path with the lowest cumulative
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
88 cost. Bowtie tries to minimize the memory impact of the descriptors, but
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
89 they can still grow very large in some cases. If you receive an error message
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
90 saying that chunk memory has been exhausted in --best mode, try adjusting
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
91 this parameter up to dedicate more memory to the descriptors. Default: 512.
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
92 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
93 parser.add_argument( '--chunkmbs', type=int, default=512 )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
94
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
95 args = parser.parse_args()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
96
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
97 # Create bismark index if necessary.
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
98 index_dir = ""
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
99 if args.own_file:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
100 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
101 Create a temporary index with the offered files from the user.
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
102 Utilizing the script: bismark_genome_preparation
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
103 bismark_genome_preparation --bowtie2 hg19/
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
104 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
105 tmp_index_dir = tempfile.mkdtemp()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
106 index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( args.own_file )[1].split( '.' )[:-1] ) )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
107 try:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
108 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
109 Create a hard link pointing to args.own_file named 'index_path'.fa.
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
110 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
111 os.symlink( args.own_file, index_path + '.fa' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
112 except Exception, e:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
113 if os.path.exists( tmp_index_dir ):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
114 shutil.rmtree( tmp_index_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
115 stop_err( 'Error in linking the reference database.\n' + str( e ) )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
116 # bismark_genome_preparation needs the complete path to the folder in which the database is stored
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
117 if args.bowtie2:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
118 cmd_index = 'bismark_genome_preparation --bowtie2 %s ' % ( tmp_index_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
119 else:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
120 cmd_index = 'bismark_genome_preparation %s ' % ( tmp_index_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
121 if args.bismark_path:
4
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
122 if os.path.exists(args.bismark_path):
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
123 # add the path to the bismark perl scripts, that is needed for galaxy
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
124 cmd_index = os.path.join(args.bismark_path, cmd_index)
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
125 else:
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
126 # assume the same directory as that script
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
127 cmd_index = 'perl %s' % os.path.join(os.path.realpath(os.path.dirname(__file__)), cmd_index)
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
128 try:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
129 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
130 tmp_stderr = open( tmp, 'wb' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
131 proc = subprocess.Popen( args=cmd_index, shell=True, cwd=tmp_index_dir, stdout=open(os.devnull, 'wb'), stderr=tmp_stderr.fileno() )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
132 returncode = proc.wait()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
133 tmp_stderr.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
134 # get stderr, allowing for case where it's very large
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
135 tmp_stderr = open( tmp, 'rb' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
136 stderr = ''
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
137 buffsize = 1048576
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
138 try:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
139 while True:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
140 stderr += tmp_stderr.read( buffsize )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
141 if not stderr or len( stderr ) % buffsize != 0:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
142 break
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
143 except OverflowError:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
144 pass
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
145 tmp_stderr.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
146 if returncode != 0:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
147 raise Exception, stderr
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
148 except Exception, e:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
149 if os.path.exists( tmp_index_dir ):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
150 shutil.rmtree( tmp_index_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
151 stop_err( 'Error indexing reference sequence\n' + str( e ) )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
152 index_dir = tmp_index_dir
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
153 else:
4
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
154 # bowtie path is the path to the index directory and the first path of the index file name
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
155 index_dir = os.path.dirname( args.index_path )
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
156
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
157 # Build bismark command
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
158 tmp_bismark_dir = tempfile.mkdtemp()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
159 output_dir = os.path.join( tmp_bismark_dir, 'results')
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
160 cmd = 'bismark %(args)s --temp_dir %(tmp_bismark_dir)s -o %(output_dir)s --quiet %(genome_folder)s %(reads)s'
4
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
161
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
162 if args.fasta:
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
163 # he query input files (specified as mate1,mate2 or singles) are FastA
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
164 cmd = '%s %s' % (cmd, '--fasta')
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
165 elif args.fastq:
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
166 cmd = '%s %s' % (cmd, '--fastq')
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
167
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
168 if args.bismark_path:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
169 # add the path to the bismark perl scripts, that is needed for galaxy
4
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
170 if os.path.exists(args.bismark_path):
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
171 cmd = os.path.join(args.bismark_path, cmd)
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
172 else:
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
173 # assume the same directory as that script
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
174 cmd = 'perl %s' % os.path.join(os.path.realpath(os.path.dirname(__file__)), cmd)
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
175
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
176 arguments = {
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
177 'genome_folder': index_dir,
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
178 'args': '',
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
179 'tmp_bismark_dir': tmp_bismark_dir,
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
180 'output_dir': output_dir,
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
181 }
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
182
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
183 additional_opts = ''
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
184 # Set up the reads
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
185 if args.mate_paired:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
186 # paired-end reads library
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
187 reads = '-1 %s ' % ( args.mate1 )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
188 reads += ' -2 %s ' % ( args.mate2 )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
189 additional_opts += ' -I %s -X %s ' % (args.min_insert, args.max_insert)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
190 else:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
191 # single paired reads library
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
192 reads = ' %s ' % ( args.single_paired )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
193
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
194
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
195 if not args.bowtie2:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
196 # use bowtie specific options
4
427fb56f2e41 - new options
bjoern-gruening
parents: 0
diff changeset
197 #additional_opts += ' --best ' # bug in bismark, --best is not available as option. Only --non-best, best-mode is activated by default
0
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
198 if args.seed_mismatches:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
199 # --seedmms
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
200 additional_opts += ' -n %s ' % args.seed_mismatches
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
201 if args.seed_len:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
202 # --seedlen
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
203 additional_opts += ' -l %s ' % args.seed_len
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
204
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
205 # alignment options
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
206 if args.bowtie2:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
207 additional_opts += ' -p %s --bowtie2 ' % args.num_threads
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
208 if args.seed_mismatches:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
209 additional_opts += ' -N %s ' % args.seed_mismatches
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
210 if args.seed_len:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
211 additional_opts += ' -L %s ' % args.seed_len
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
212 if args.seed_extention_attempts:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
213 additional_opts += ' -D %s ' % args.seed_extention_attempts
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
214 if args.max_reseed:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
215 additional_opts += ' -R %s ' % args.max_reseed
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
216 if args.no_discordant:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
217 additional_opts += ' --no-discordant '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
218 if args.no_mixed:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
219 additional_opts += ' --no-mixed '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
220 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
221 if args.maqerr:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
222 additional_opts += ' --maqerr %s ' % args.maqerr
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
223 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
224 if args.skip_reads:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
225 additional_opts += ' --skip %s ' % args.skip_reads
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
226 if args.qupto:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
227 additional_opts += ' --qupto %s ' % args.qupto
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
228 if args.phred64:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
229 additional_opts += ' --phred64-quals '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
230 if args.suppress_header:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
231 additional_opts += ' --sam-no-hd '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
232 if args.output_unmapped_reads or ( args.output_unmapped_reads_l and args.output_unmapped_reads_r):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
233 additional_opts += ' --un '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
234 if args.output_suppressed_reads or ( args.output_suppressed_reads_l and args.output_suppressed_reads_r):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
235 additional_opts += ' --ambiguous '
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
236
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
237 arguments.update( {'args': additional_opts, 'reads': reads} )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
238
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
239 # Final command:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
240 cmd = cmd % arguments
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
241
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
242 # Run
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
243 try:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
244 tmp_out = tempfile.NamedTemporaryFile().name
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
245 tmp_stdout = open( tmp_out, 'wb' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
246 tmp_err = tempfile.NamedTemporaryFile().name
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
247 tmp_stderr = open( tmp_err, 'wb' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
248 proc = subprocess.Popen( args=cmd, shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
249 returncode = proc.wait()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
250 tmp_stderr.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
251 # get stderr, allowing for case where it's very large
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
252 tmp_stderr = open( tmp_err, 'rb' )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
253 stderr = ''
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
254 buffsize = 1048576
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
255 try:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
256 while True:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
257 stderr += tmp_stderr.read( buffsize )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
258 if not stderr or len( stderr ) % buffsize != 0:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
259 break
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
260 except OverflowError:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
261 pass
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
262 tmp_stdout.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
263 tmp_stderr.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
264 if returncode != 0:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
265 raise Exception, stderr
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
266
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
267 # TODO: look for errors in program output.
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
268 except Exception, e:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
269 stop_err( 'Error in bismark:\n' + str( e ) )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
270
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
271
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
272 # collect and copy output files
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
273 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
274 if args.output_report_file:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
275 output_report_file = open(args.output_report_file, 'w+')
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
276 for line in fileinput.input(glob( os.path.join( output_dir, '*.txt') )):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
277 output_report_file.write(line)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
278 output_report_file.close()
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
279 """
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
280
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
281 if args.output_suppressed_reads:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
282 shutil.move( glob(os.path.join( output_dir, '*ambiguous_reads.txt'))[0], args.output_suppressed_reads )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
283 if args.output_suppressed_reads_l:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
284 shutil.move( glob(os.path.join( output_dir, '*ambiguous_reads_1.txt'))[0], args.output_suppressed_reads_l )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
285 if args.output_suppressed_reads_r:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
286 shutil.move( glob(os.path.join( output_dir, '*ambiguous_reads_2.txt'))[0], args.output_suppressed_reads_r )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
287
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
288 if args.output_unmapped_reads:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
289 shutil.move( glob(os.path.join( output_dir, '*unmapped_reads.txt'))[0], args.output_unmapped_reads )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
290 if args.output_unmapped_reads_l:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
291 shutil.move( glob(os.path.join( output_dir, '*unmapped_reads_1.txt'))[0], args.output_unmapped_reads_l )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
292 if args.output_unmapped_reads_r:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
293 shutil.move( glob(os.path.join( output_dir, '*unmapped_reads_2.txt'))[0], args.output_unmapped_reads_r )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
294
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
295 shutil.move( glob( os.path.join( output_dir, '*.sam'))[0] , args.output)
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
296
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
297 # Clean up temp dirs
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
298 if args.own_file:
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
299 if os.path.exists( tmp_index_dir ):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
300 shutil.rmtree( tmp_index_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
301 if os.path.exists( tmp_bismark_dir ):
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
302 shutil.rmtree( tmp_bismark_dir )
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
303
36d124f44c0a inital commit
bjoern-gruening
parents:
diff changeset
304 if __name__=="__main__": __main__()