Discard reads from file that fail to map as a proper pair.

[paolo.kunder]

Dear All,
I would like use Samtools to discard reads from the bam file that fail to map as a proper pair.

Is this command correct?

samtools view -b-h -f 4 -o out.bam in.bam

I am wondering about the integer after the f option

Thanks,
Paolo

[syfo]

Dear All,
I would like use Samtools to discard reads from the bam file that fail to map as a proper pair.

Is this command correct?

samtools view -b-h -f 4 -o out.bam in.bam

I am wondering about the integer after the f option

Thanks,
Paolo

I would have done:

samtools view -bh -f2 in.bam > out.bam

-f4 gets you the unmapped reads
FYI, a flag reader is available here

[swbarnes2]

You want

samtools view -bhF 2 in.bam > out.bam

You want to filter out all the reads that didn't map right, so the F needs to be uppercase.

[syfo]

You want

samtools view -bhF 2 in.bam > out.bam

You want to filter out all the reads that didn't map right, so the F needs to be uppercase.

mhm.. isn't it precisely the other way around? I think that "-f2" gives you proper pairs and "-F2" filters them out.

[winsettz]

From samtools (http://samtools.sourceforge.net/samtools.shtml)

-f INT Only output alignments with all bits in INT present in the FLAG field. INT can be in hex in the format of /^0x[0-9A-F]+/ [0]
-F INT Skip alignments with bits present in INT [0]

-f to output only a particular bit, -F to skip.

From the SAM spec: (http://samtools.sourceforge.net/SAM1.pdf)

FLAG: bitwise FLAG. Each bit is explained in the following table:
Bit Description
0x1 template having multiple segments in sequencing
0x2 each segment properly aligned according to the aligner
0x4 segment unmapped

SAM spec is not as clear as I would like. 4 is obvious, and I've thought of 2 as meaning "mapped to one position" and 1 meaning "maps in multiple regions".

Though the flag explainer is new to me; thanks for posting it.

[swbarnes2]

You want

samtools view -bhF 2 in.bam > out.bam

You want to filter out all the reads that didn't map right, so the F needs to be uppercase.

Edit: Oops, brain fart, I was wrong.

[syfo]

Even better: add -F12 to make sure you remove the ones that are not mapped and whose mate is not mapped. Believe it or not, the "mapped in proper pair" flag can be given to... unmapped reads.

Explanation from the SAM format specs (according to a previous post): "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.".

Conclusion: don't trust the sole "proper pair" nor the "mapq" flag and always control for the mapping status.