Skip to content

Commit

Permalink
Fix two issues in using bwa mem for alignments (#633)
Browse files Browse the repository at this point in the history
* Fix bug in setting bwa mem options in --sensitive mode

* Remove alignment score filter on bwa mem

Based on bwa's definition of its alignment score, it
might be reasonable for it to output an alignment score
less than the specified `-T` value.
  • Loading branch information
haydenm authored and tomkinsc committed Jun 9, 2017
1 parent 93b7061 commit dd91a30
Show file tree
Hide file tree
Showing 2 changed files with 4 additions and 15 deletions.
2 changes: 1 addition & 1 deletion reports.py
Original file line number Diff line number Diff line change
Expand Up @@ -746,7 +746,7 @@ def align_and_plot_coverage(

bwa_opts = aligner_options.split()
if sensitive:
bwa_opts + "-k 12 -A 1 -B 1 -O 1 -E 1".split()
bwa_opts += "-k 12 -A 1 -B 1 -O 1 -E 1".split()

# get the quality threshold from the opts
# for downstream filtering
Expand Down
17 changes: 3 additions & 14 deletions tools/bwa.py
Original file line number Diff line number Diff line change
Expand Up @@ -156,29 +156,18 @@ def align_mem_one_rg(self, inBam, refDb, outBam, rgid=None, options=None, min_qu
readgroup_line = line

assert len(readgroup_line) > 0
aln_bam_prefilter = util.file.mkstempfname('.prefiltered.bam')

tmp_bam_aligned = util.file.mkstempfname('.aligned.bam')
# rather than reheader the alignment bam file later so it has the readgroup information
# from the original bam file, we'll pass the RG line to bwa to write out
self.mem(one_rg_inBam, refDb, aln_bam_prefilter, options=options+['-R', readgroup_line.rstrip("\n").rstrip("\r")], min_qual=min_qual, threads=threads)
self.mem(one_rg_inBam, refDb, tmp_bam_aligned, options=options+['-R', readgroup_line.rstrip("\n").rstrip("\r")], min_qual=min_qual, threads=threads)

# if there was more than one RG in the input, we had to create a temporary file with the one RG specified
# and we can safely delete it this file
# if there was only one RG in the input, we used it directly and should not delete it
if removeInput:
os.unlink(one_rg_inBam)

# @haydenm says:
# For some reason (particularly when the --sensitive option is on), bwa
# doesn't listen to its '-T' flag and outputs alignments with score less
# than the '-T 30' threshold. So filter these:
if min_qual > 0:
tmp_bam_aligned = util.file.mkstempfname('.aligned.bam')
tools.samtools.SamtoolsTool().view(["-b", "-h", "-q", str(min_qual)], aln_bam_prefilter, tmp_bam_aligned)
os.unlink(aln_bam_prefilter)
else:
shutil.move(aln_bam_prefilter, tmp_bam_aligned)

# if the aligned bam file contains no reads after filtering
# just create an empty file
if tools.samtools.SamtoolsTool().count(tmp_bam_aligned) == 0:
Expand Down

0 comments on commit dd91a30

Please sign in to comment.