This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BWA missing reads in output

I'm having trouble getting BWA mem (v0.7.12) to consistently output all reads in the input fastqs. I'm working on a simulation study, and I'm generating smallish sets of fake reads, then aligning them with bwa. The paired input fastqs have a few thousand reads each. The missing reads are all from the ends of the input fastq files as far as I can tell, and subsetting the input fastqs to focus on reads missing from a previous attempt can make the reads show back up again. It looks a little like under some circumstances BWA is forgetting to flush a buffer, or something like that.

Since I'm simulating reads I know exactly where they should go, and I give them each unique ids so I can verify which ones exist in the fastqs and the output sams. The behavior seems identical on both linux and mac.

Are there circumstances under which BWA will refuse to align / output a read, or set of reads? Everything in the fastqs should end up in the .sam, right? Anyone else notice anything similar?

bwa alignment
• 3,479 views • link
updated by Ram Ram • written by brendanofallon • brendanofallon

output wc -l of fastq and samtools flagstat of bam file.

• 0 views • link
written by apelin20 apelin20

Here's a typical run for me:

bofallon@limital:/tmp-working/$ grep '^@' input_r1.fq | wc -l
1250

bofallon@limital:/tmp-working/$ grep '^@' input_r2.fq | wc -l
1250

bofallon@limital:/tmp-working/$ samtools flagstat input_r1.fq.bam
  1222 + 0 in total (QC-passed reads + QC-failed reads)
  0 + 0 secondary
  0 + 0 supplementary
  0 + 0 duplicates
  1222 + 0 mapped (100.00%:-nan%)
  1222 + 0 paired in sequencing
  611 + 0 read1
  611 + 0 read2
  1222 + 0 properly paired (100.00%:-nan%)
  1222 + 0 with itself and mate mapped
  0 + 0 singletons (0.00%:-nan%)
  0 + 0 with mate mapped to a different chr
  0 + 0 with mate mapped to a different chr (mapQ>=5)

1250 input read pairs... and 1222 output read pairs.

(I should mention that with this simulated data no base qualities are '@', so the grepping above accurately counts reads - this won't be true for real data)

• 0 views • link
updated by Ram Ram • written by brendanofallon • brendanofallon

What command line(s) did you use to go from fastq to bam?

• 0 views • link
written by Sean Davis Sean Davis

0 answers

No answers yet.

Log in to answer this question.

More posts like this