This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Separating unmapped reads from bam files is not working using samtools -f 4

I mapped reads using bwa mem and created a sorted bam file. An example flagstat is as follows: 

12446425 + 0 in total (QC-passed reads + QC-failed reads)
178116 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
10794644 + 0 mapped (86.73% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

So mapping is at 86%. However, when I try to separate out the unmapped reads using this code: 

samtools view -f 4 -H -b in.bam > unmapped.bam

I get no reads. 

0 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 mapped (N/A : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

What am I missing?

genome samtools bwa
• 1,662 views • link
updated by GenoMax GenoMax • written by sorrymouse sorrymouse

1 answer 

  -H       print SAM header only (no alignments)

Leave out the -H, you probably mean -h to print a file with header but -b implies -h therefore:

samtools view -f 4 -b in.bam > unmapped.bam
• 0 views • link
written by ATpoint ATpoint

I seem to have had more than one problem here - YES this would very much be an issue. I can't believe I missed that.

Also, not including the -H is still producing no unmapped reads.

For some reason, going back to the sam file, rather than the bam, is working. No idea why.

Thanks for noticing that coding error, clearly nothing was going to work as long as that mistake was there.

• 0 views • link
written by sorrymouse sorrymouse

Thanks for noticing that coding error, clearly nothing was going to work as long as that mistake was there.

Don't worry, these little mistaken happen to all of us at times :)

• 0 views • link
written by ATpoint ATpoint

Log in to answer this question.

More posts like this