No announcement yet.

BWA output

  • Filter
  • Time
  • Show
Clear All
new posts

  • BWA output

    Hi all,

    I'm trying to screen paired Illumina RNAseq reads against a database of putative contaminants prior to assembly. BWA (version 7.0.3a, in case it matters) ran beautifully. Now I'm trying to extract the names of all the reads that mapped to something in my contaminant database. Basically, if a read OR its mate matched a putative contaminant, I want to toss it out. A colleague told me that the output is easy to parse, given that all the unmapped reads will have a "*" in the third column of the output line. I wrote up a quick perl script to parse the output based on the "*".

    Later, I learned that its possible to convert BWA's sam output to bam (via samtools) and then count the mapped reads with the following command:
    samtools view -c -F4 <bam file>
    (I found that here:

    The issue is that I get two different answers, and I can't figure out why. What's wrong with my approach?

  • #2
    I re-read this and wonder if the description of my perl script is unclear. I basically counted anything that did NOT have an "*" in the 3rd column. I just didn't want anyone to think I was getting inconsistent results because I had this backwards.

    Thanks in advance for your feedback!