-
Notifications
You must be signed in to change notification settings - Fork 2
/
align.py
executable file
·369 lines (322 loc) · 12.5 KB
/
align.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
#!/usr/bin/env python3
import argparse
import distutils.spawn
from distutils.version import LooseVersion
import logging
import os
import pathlib
import subprocess
import sys
assert sys.version_info.major >= 3, 'Python 3 required'
ALIGNER_DATA = {
'bowtie2': {
'required': ('bowtie2-build', 'bowtie2'),
'index-endings': ('1.bt2', '2.bt2', '3.bt2', '4.bt2', 'rev.1.bt2', 'rev.2.bt2'),
'opts': [],
},
'bwa': {
'required': ('bwa',),
'index-endings': ('amb', 'ann', 'bwt', 'pac', 'sa'),
'opts': ['-M'],
},
}
DESCRIPTION = """Align reads.
This will automatically run all the commands needed to get you from raw reads to an indexed BAM
file, including indexing the reference sequence."""
def make_argparser():
parser = argparse.ArgumentParser(description=DESCRIPTION)
parser.add_argument('aligner', choices=ALIGNER_DATA.keys(),
help="Aligner to use. 'bwa' will use BWA-MEM.")
parser.add_argument('ref', metavar='ref.fa', type=pathlib.Path,
help='Reference sequence.')
parser.add_argument('reads1', metavar='reads_1.fq', type=pathlib.Path,
help='Reads, mate 1.')
parser.add_argument('reads2', metavar='reads_2.fq', type=pathlib.Path,
help='Reads, mate 2.')
parser.add_argument('-o', '--out', type=pathlib.Path,
help='Output file. The aligned BAM will be written to this path.')
parser.add_argument('-f', '--format', choices=('sam', 'bam'),
help="Output format. If this option isn't given, and a filename ending in '.sam' is given to "
'--out, it will output SAM. Otherwise it will default to BAM.')
parser.add_argument('-c', '--clobber', action='store_true',
help='Overwrite intermediate and output files without prompting. Otherwise, the script will '
'fail if one of these files already exists.')
parser.add_argument('-R', '--ref-base',
help='The base path for the reference index. E.g. "-R path/to/ref" if your index files are '
'named like "path/to/ref.bwt" and "path/to/ref.pac". Default: the same path as the actual '
'reference file.')
parser.add_argument('-i', '--keep-index', action='store_true',
help="Don't delete the reference index files. Default: leave it as you found it "
"(clean up if they didn't exist, leave them if they did).")
parser.add_argument('-I', '--delete-index', action='store_true',
help='Delete the reference index files, even if they already existed.')
parser.add_argument('-t', '--threads', type=int, default=1,
help='Threads to use when aligning. For bowtie2, this also speeds up indexing. For modern '
'versions of samtools (>= 1.0), this speeds up sorting. Default: %(default)s')
opts_str = ''
for aligner, data in ALIGNER_DATA.items():
opts = data.get('opts')
if opts:
if opts_str:
opts_str += ', '
opts_str += "'"+' '.join(opts)+f"' for '{aligner}'"
parser.add_argument('-O', '--aligner-opts', type=split_opt_list,
help="Options to pass to the alignment command. If giving a single option, prepend '|' to "
"keep it from being parsed as an argument to this script. E.g. \"-O '|-M'\". "
'Defaults: '+opts_str+', none for the rest.')
parser.add_argument('-N', '--name-sort', dest='sort_key', action='store_const', const='name',
default='coord',
help='Sort the output BAM by name instead of coordinate.')
parser.add_argument('-l', '--log', type=argparse.FileType('w'), default=sys.stderr,
help='Print log messages to this file instead of to stderr. Warning: Will overwrite the file.')
volume = parser.add_mutually_exclusive_group()
volume.add_argument('-q', '--quiet', dest='volume', action='store_const', const=logging.CRITICAL,
default=logging.WARNING)
volume.add_argument('-v', '--verbose', dest='volume', action='store_const', const=logging.INFO)
volume.add_argument('-D', '--debug', dest='volume', action='store_const', const=logging.DEBUG)
return parser
def main(argv):
parser = make_argparser()
args = parser.parse_args(argv[1:])
logging.basicConfig(stream=args.log, level=args.volume, format='%(message)s')
# Check for required commands.
missing = []
for command in ALIGNER_DATA[args.aligner]['required']:
if not distutils.spawn.find_executable(command):
missing.append(command)
if missing:
fail("Error: Missing required command(s) '"+"', '".join(missing)+"'.")
# Determine paths.
if args.ref_base:
ref_base = args.ref_base
else:
ref_base = str(args.ref)
format = get_format(args.out, args.format)
sam_path, out_path = get_paths(args.reads1, args.out, format)
if not args.clobber:
for path in sam_path, out_path:
if path.exists():
fail(f'Error: output path {path} already exists. Use -c/--clobber to overwrite.')
was_indexed = is_indexed(args.aligner, ref_base)
# Index reference (if needed).
if not was_indexed:
clear_indices(args.aligner, ref_base)
index_ref(args.aligner, args.ref, ref_base, threads=args.threads)
# Do alignment.
align(
args.aligner, ref_base, args.reads1, args.reads2, sam_path, threads=args.threads,
opts=args.aligner_opts
)
if format == 'bam':
# Convert output.
convert(sam_path, out_path, sort_key=args.sort_key, threads=args.threads)
# Index output.
if args.sort_key == 'coord':
index_bam(out_path)
# Clean up intermediate file(s).
os.remove(sam_path)
if (not was_indexed and not args.keep_index) or args.delete_index:
clear_indices(args.aligner, ref_base)
if out_path.is_file():
logging.error(f'Success! Output is in {str(out_path)!r}')
else:
fail(f'Error: Output file missing {str(out_path)!r}')
def get_paths(reads1_arg, out_arg, format):
base = get_reads_base(reads1_arg)
if format == 'sam':
if out_arg:
sam_path = out_arg
else:
sam_path = pathlib.Path(base+'.sam')
out_path = sam_path
elif format == 'bam':
sam_path = pathlib.Path(base+'.sam')
if out_arg:
out_path = out_arg
else:
out_path = pathlib.Path(base+'.bam')
return sam_path, out_path
def get_reads_base(reads_path: pathlib.Path) -> str:
"""Get the basename of a file, omitting any `_1` or `_2` before the extension.
Takes a path, possibly including directories, and returns a string of the whole
path, but with the file extension and any `_1` or `_2` removed."""
basename = reads_path.stem
if basename.endswith('_1') or basename.endswith('_2'):
basename = basename[:-2]
return str(reads_path.parent / basename)
def get_format(out_arg, format_arg):
if format_arg:
return format_arg
if out_arg is not None:
ext = out_arg.suffix.lower()
if ext == '.bam':
return 'bam'
elif ext == '.sam':
return 'sam'
# Default: BAM
return 'bam'
def is_indexed(aligner, ref_base):
for ending in ALIGNER_DATA[aligner]['index-endings']:
ref_path_str = ref_base+'.'+ending
if os.path.exists(ref_path_str):
if not os.path.isfile(ref_path_str):
fail(f'Error: Index path {ref_path_str!r} exists, but is not a regular file.')
else:
return False
return True
def clear_indices(aligner, ref_base):
for ending in ALIGNER_DATA[aligner]['index-endings']:
ref_path_str = ref_base+'.'+ending
if os.path.isfile(ref_path_str):
os.remove(ref_path_str)
def index_ref(aligner, ref, ref_base, threads=1):
kwargs = {}
if aligner == 'bowtie2':
# $ bowtie2-build ref.fa ref
cmd = ['bowtie2-build']
version = get_bowtie2_version()
if version is None:
logging.warning('Warning: Unable to determine bowtie2 version.')
elif version >= LooseVersion('2.2.7'):
cmd.extend(['--threads', threads])
cmd.extend([ref, ref_base])
kwargs['stdout'] = subprocess.DEVNULL
elif aligner == 'bwa':
# $ bwa index -a algo -p ref ref.fa
# Use the 'is' indexing algorithm, unless the reference is > 2GB.
if os.path.getsize(ref) < 2000000000:
algorithm = 'is'
else:
algorithm = 'bwtsw'
cmd = ('bwa', 'index', '-a', algorithm, '-p', ref_base, ref)
run_command(cmd, **kwargs)
def align(aligner, ref_base, reads1_path, reads2_path, sam_path, threads=1, opts=None):
if opts is None:
opts = ALIGNER_DATA[aligner]['opts']
if aligner == 'bowtie2':
# $ bowtie2 -x ref -1 reads_1.fq -2 reads_2.fq -S align.sam
cmd = (
['bowtie2'] + opts +
['-p', threads, '-x', ref_base, '-1', reads1_path, '-2', reads2_path, '-S', sam_path]
)
run_command(cmd)
elif aligner == 'bwa':
# $ bwa mem [opts] ref reads_1.fq reads_2.fq > align.sam
cmd = ['bwa', 'mem'] + opts + ['-t', threads, ref_base, reads1_path, reads2_path]
with sam_path.open('wb') as sam_file:
run_command(cmd, stdout=sam_file)
def convert(sam_path, bam_path, sort_key='coord', threads=1):
# $ samtools view -Sb align.sam | samtools sort -o - dummy > align.bam
# --- or ---
# $ samtools view -Sb align.sam | samtools sort -O bam - > align.bam
version = get_samtools_version()
if version is None:
logging.warning('Warning: Unable to determine samtools version.')
cmd1 = ['samtools', 'view', '-Sb', sam_path]
# New samtools sort call format introduced in 0.2.0-rc9.
if version is not None and version < LooseVersion('0.2.0-rc9'):
cmd2 = ['samtools', 'sort', '-o', '-', 'dummy']
else:
cmd2 = ['samtools', 'sort', '-O', 'bam', '-']
# -@ argument introduced sometime between 0.1.19 and 1.0.
if version is not None and version > LooseVersion('1.0'):
cmd2[2:2] = ['-@', threads]
if sort_key == 'name':
cmd2[2:2] = ['-n']
cmd1 = [str(arg) for arg in cmd1]
cmd2 = [str(arg) for arg in cmd2]
logging.warning(
'$ '+' '.join(cmd1)+' \\\n'
+' | '+' '.join(cmd2)+' \\\n'
+' > '+str(bam_path)
)
proc1 = subprocess.Popen(cmd1, stdout=subprocess.PIPE)
with bam_path.open('wb') as bam_file:
proc2 = subprocess.Popen(cmd2, stdin=proc1.stdout, stdout=bam_file)
proc1.stdout.close()
for cmd, proc in zip((cmd1, cmd2), (proc1, proc2)):
if proc.wait():
fail(f'Error: {cmd[0]} {cmd[1]} returned with a non-zero exit code.')
def index_bam(bam_path):
# $ samtools index align.bam
index_path = bam_path.parent / (bam_path.name+'.bai')
if index_path.exists():
os.remove(index_path)
cmd = ('samtools', 'index', bam_path)
if bam_path.is_file():
run_command(cmd)
else:
fail(f'Error: Output BAM missing ({bam_path})')
def get_samtools_version(exe='samtools'):
cmd = (exe,)
result = subprocess.run(cmd, stdout=subprocess.DEVNULL, stderr=subprocess.PIPE)
if result.returncode != 1:
return None
output = str(result.stderr, 'utf8')
# Find the version line.
line_fields = None
for line in output.splitlines():
if line.lower().startswith('version:'):
line_fields = line.split()
if line_fields is None:
return None
# Find the version number in the line.
ver_str = line_fields[1]
# Verify it looks like a version number: does it start with two decimal-separated integers?
ver_fields = ver_str.split('.')
if len(ver_fields) <= 1:
return None
try:
int(ver_fields[0])
int(ver_fields[1])
except ValueError:
return None
logging.info(f'Info: Successfully determined samtools version to be {ver_str}.')
return LooseVersion(ver_str)
def get_bowtie2_version(exe='bowtie2-build'):
cmd = (exe, '--version')
result = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL)
if result.returncode != 0:
return None
output = str(result.stdout, 'utf8')
# Find the version line.
line_fields = None
for line_num, line in enumerate(output.splitlines()):
if line_num == 0:
line_fields = line.split()
if not line_fields:
continue
exe_path = pathlib.Path(line_fields[0])
if exe_path.name.startswith('bowtie2-build'):
break
if line_fields is None:
return None
# Find the version number in the line.
ver_str = line_fields[-1]
# Verify it looks like the right version number:
# Is it at least two decimal-delimited fields, the first of which is 2?
ver_fields = ver_str.split('.')
if len(ver_fields) > 1 and ver_fields[0] == '2':
logging.info(f'Info: Successfully determined bowtie2 version to be {ver_str}.')
return LooseVersion(ver_str)
else:
return None
def run_command(cmd, **kwargs):
cmd_strs = [str(arg) for arg in cmd]
logging.warning('$ '+' '.join(cmd_strs))
if subprocess.call(cmd_strs, **kwargs):
fail(f'Error: {cmd[0]} returned with a non-zero exit code.')
def split_opt_list(opts_str):
opts_str2 = opts_str.lstrip('|')
return opts_str2.split()
def fail(message):
logging.critical(message)
if __name__ == '__main__':
sys.exit(1)
else:
raise Exception('Unrecoverable error')
if __name__ == '__main__':
try:
sys.exit(main(sys.argv))
except BrokenPipeError:
pass