Separating unmapped reads from bam files is not working using samtools -f 4

Question

[sorrymouse]

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?

Answer 1

[ATpoint]

  -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

[sorrymouse]

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.

[ATpoint]

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 :)