annotate bismark_wrapper/bismark_wrapper.py @ 1:183de9d00131 draft

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