Strange Behavior Filtering Reads With Samtools

Question

[Linda]

I am trying to filter the results of a paired end alignment with samtools. I would like to get all reads where both the pairs are unaligned. I tried

samtools view  -f 12 R28.bam | cut -f1-9|head
HWUSI-EAS570R_0028:2:96:4625:6148#0    141    Contig18298    1907    0    100M    =    1907    0
HWUSI-EAS570R_0028:1:86:7595:20840#0    141    Contig18298    1916    25    100M    =    1916    0
HWUSI-EAS570R_0028:3:120:15625:9857#0    93    Contig18301    2004    25    100M    =    2004    0
HWUSI-EAS570R_0028:5:89:9474:1987#0    157    Contig18321    1910    37    6M2D2M1I91M    =    1910    0
HWUSI-EAS570R_0028:5:53:6864:4533#0    141    Contig18345    1921    25    100M    =    1921    0
HWUSI-EAS570R_0028:1:27:10901:2866#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:1:36:7563:3669#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:1:37:12448:8115#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:2:68:9495:3274#0    157    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:3:20:9943:7308#0    93    Contig18367    1918    0    100M    =    1918    0

The 100M indicates to me that these are matching the reference sequence. So I tried a two step filter.

samtools view -b -f 4 R28.bam | samtools view -f 8 - |cut -f1-9|head

HWUSI-EAS570R_0028:2:96:4625:6148#0    141    Contig18298    1907    0    100M    =    1907    0
HWUSI-EAS570R_0028:1:86:7595:20840#0    141    Contig18298    1916    25    100M    =    1916    0
HWUSI-EAS570R_0028:3:120:15625:9857#0    93    Contig18301    2004    25    100M    =    2004    0
HWUSI-EAS570R_0028:5:89:9474:1987#0    157    Contig18321    1910    37    6M2D2M1I91M    =    1910    0
HWUSI-EAS570R_0028:5:53:6864:4533#0    141    Contig18345    1921    25    100M    =    1921    0
HWUSI-EAS570R_0028:1:27:10901:2866#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:1:36:7563:3669#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:1:37:12448:8115#0    93    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:2:68:9495:3274#0    157    Contig18367    1918    0    100M    =    1918    0
HWUSI-EAS570R_0028:3:20:9943:7308#0    93    Contig18367    1918    0    100M    =    1918    0

Now, this is strange, because I thought the first filter should have removed all reads that aligned anywhere? So what is going on here? My aim is to collect all unaligned pairs.

Answer 1

[Istvan Albert]

The SAM specification states that:

Bit 0x4 is the only reliable place to tell whether the segment is unmapped. If 0x4 is set, no assumptions can be made about RNAME, POS, CIGAR, MAPQ, bits 0x2, 0x10 and 0x100 and the bit 0x20 of the next segment in the template.

This means that you don't need to worry about the other fields since those are not meaningful. These are probably not cleared due to some optimization reason.